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NONLINEAR PROBABILISTIC FINITE ELEMENT MODELS OF 
LAMINATED COMPOSITE SHELLS 

S. P. Engelstad and J. N. Reddy 

Department of Engineering Science and Mechanics 
Virginia Polytechnic Institute and State University, Blacksburg, VA 24061 

SUMMARY 

A probabilistic finite element analysis procedure for laminated composite 
shells has been developed. A total Lagrangian finite element formulation, 
employing a degenerated 3— D laminated composite shell element with the full 
Green— Lagrange strains and first— order shear deformable kinematics, forms the 
modeling foundation. The first-order second-moment technique for 
probabilistic finite element analysis of random fields is employed and results are 
presented in the form of mean and variance of the structural response. The 
effects of material nonlinearity are included through the use of a 
rate— independent anisotropic plasticity formulation with the macroscopic point 
of view. Both ply-level and micromechanics— level random variables can be 
selected, the latter by means of the Aboudi micromechanics model. A number 
of sample problems are solved to verify the accuracy of the procedures 
developed and to quantify the variability of certain material type/structure 
combinations. Experimental data is compared in many cases, and the Monte 
Carlo simulation method is used to check the probabilistic results. In general, 
the procedure is quite effective in modeling the mean and variance response of 
the linear and nonlinear behavior of laminated composite shells. 
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1. INTRODUCTION 


1.1 Motivation 

The use of fiber reinforced composite materials in modern engineering 
structural design has become a common practice. Organic matrix composite 
materials such as graphite— epoxy have been used extensively with substantial 
weight savings along with additional benefits such as dimensional stability. 
More advanced materials such as metal matrix composites are being developed 
for use in aerospace structures where temperature requirements exceed the 
limits of typical organic matrix composites. Analytical methods for metal 
matrix composites are a "hot" item in current research, and in order to use 
these materials, methods for determining design limitations are necessary. 

In the area of structural analysis, the finite element method (FEM) has 
become the most widely used analysis tool. Developments and improvements 
have progressed over the last two decades to the point where many families of 
reliable finite element codes are in place. For a typical FEM analysis, it is a 
fair statement to say that uncertainties caused by modeling inaccuracies are far 
outweighed by uncertainties in the material properties, geometry, and loading of 
the problem. This is the primary reason for the development of probabilistic 
finite element methods, so that these uncertainties can be modeled within the 
framework of statistical methods. 

The benefits from the use of statistical methods in design are numerous. 
They include: 1) a reduction in design conservatism, so that materials can be 
used to their true capacity; 2) the ability to quantify the amount of variability 
itself, and 3) the estimation of the risk of exceeding design variables. As 
opposed to a deterministic yes— no failure analysis, risk of failure or "reliability" 
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of a structure can be much more meaningful to the designer. This reliability 
can only be determined if the uncertainties of the problem are included in the 
model. 

As composite materials are incorporated into modern structural 
components, lack of an experience base induces a trend towards conservatism in 
the design. Thus much of the available savings in weight and long term costs 
are lost. Since more design variables exist when composites are involved, and 
the manufacturing processes for producing composite materials themselves are 
more complex, then more variability can exist in a design produced with 
composites versus conventional materials. Thus the motivation of this research 
is to develop methods for probabilistic structural analysis of composite 
materials. Laminated shell type structures are selected as the focus, and both 
geometric and material nonlinearities are included. 

In the following section, a summary of the selected methods used in the 
present study, along with some discussions involving orginality are presented. 

1.2 Present Study 

The main objective of this study is to develop a probabilistic finite 
element procedure to be used for reliability analysis of laminated composite 
shells. Since these structures often exhibit both geometric and material 
nonlinear behavior, these nonlinearities have been included in the development. 
The discussion to follow describes the deterministic and probabilistic methods 
incorporated into the computational procedure. 
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A degenerated 3— D shell element formulation, incorporating the 
first-order shear deformable kinematics, has been chosen to model laminated 
composite shell structures. For organic matrix composites, the importance of 
shear deformation has been covered in the literature quite extensively. An 
example has been given in this work which illustrates this importance in the 
postbuckling range, both for determining postbuckling response and failure. 

The second— moment technique for probabilistic finite element analysis 
has been employed. The random variables built into the model include the ply 
stiffnesses, orientation angles, and ply thicknesses. Using a micromechanics 
theory, micromechanics level random variables such as fiber and matrix 
stiffnesses and volume ratios can be selected. 

Geometric nonlinearity is based on the total Lagrangian approach with 
the full Green— Lagrange strains. Material nonlinearity is incorporated using 
classical rate— independent plasticity and a general orthotropic yield function. 
The radial return algorithm for plane stress has been extended to calculate the 
elastic— plastic stresses with the orthotropic yield function. 

The probabilistic responses in the computer program, developed during 
this study, are the first and second probabilistic moments of the structural 
responses, which include deflection, strain, and stress. It was decided that 
verification of the second— moment methods using Monte Carlo procedures 
would be made easier if these moments (mean and variance) were calculated. 
Also, the variability of the response due to individual random variables can be 
quantified (sensitivity analysis) by estimating the variance. In order to become 
familiar with all the computational subtleties in the probabilistic finite element 
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method, the mean and variance calculations are the natural first step. 
Obviously, future work involves incorporation of the reliability estimation 
techniques. 

The originality in this study lies in the application of the 
second-moment method for probabilistic finite element analysis to geometric 
and material nonlinear composite shells. Previous work in the literature for 
composites only involved linear analysis of plates, in which classical lamination 
theory was used. Other works involved geometric or material nonlinearities in a 
probabilistic finite element format, but only for isotropic materials and usually 
only for plates. It is of interest here to analyze more realistic laminated 
composite structures, so the shell element with shear deformation theory has 
been employed. The computational difficulties in analyzing actual laminated 
composite panels, involving models with large numbers of layers and degrees of 
freedom, and also that proceed deep into the postbuckling range, have been 
investigated. Several comparisons have been made with experimental results as 
well. 

A review of the literature involving several areas of research of 
importance to this study is presented in Sec. 2. Section 3 contains the 
theoretical development of the degenerated 3-D shell element for laminated 
composite shells and the finite element formulation of the incremental equations 
of motion including geometric nonlinearity. The Aboudi micromechanics theory 
is developed in Sec. 4, and the anisotropic plasticity formulation is 
presented in Sec. 5. Section 6 contains the development of the 
second— moment probabilistic finite element method, including discussions on 
computational saving techniques. A number of illustrative problems are solved 
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in Sec . 7, which demonstrate various aspects and capabilities of the overall 
procedure. Section 8 provides a summary, conclusions, and recommendations 
for future work. 

Throughout this study it has been assumed that the reader has a 
basic understanding of probability concepts and terms. However, to aid in this 
area, a review of reliability estimation theory, along with a brief explanation of 
basic terms is presented in Appendix A. Appendix B describes the Monte Carlo 
simulation technique. 
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2. LITERATURE REVIEW 

The purpose of this review is to assess the current state of probabilistic 
finite element methods in general, and in particular to include any work in 
relation to composite materials. A brief discussion of various micromechanics 
models, macroscopic anisotropic plasticity theories, and reliability estimation 
methods will be given as well. 

2.1 Probabilistic Finite Element Methods 

Probabilistic finite element methods have generally evolved along two 
major paths [33]. The first one uses standard statistical methods in conjunction 
with the finite element code. The second method involves non— statistical 
techniques and differs little from deterministic methods. These two categories 
will be discussed, with emphasis on the second one since it contains the primary 
concentration of effort in current research. 

2.1.1 Statistical Approach 

Typically these methods involve statistical sampling procedures known 
as Monte Carlo simulation. Sampling from a known multivariate distribution 
function is conducted and due to the ’weak law of large numbers’ [1], large 
sample sizes must be used in order to converge to the approximately correct 
statistical parameters of the response. These methods therefore become very 
expensive, as finite element solutions must be produced for each sample. When 
correlation of the random variables exist, transformations must be performed 
prior to simulation since simulation involves independent sampling. In order to 
cut down on the number of samples required for the direct Monte Carlo 


6 


simulation technique, other methods of sampling such as stratified sampling and 
Latin hypercube sampling are often used [1—3]. Since the literature in this 
category is quite extensive, a few of the most important and widely used 
statistical methods are reviewed. 

An example of the use of Monte Carlo techniques to study a spatial 
stochastic process is given by Ma and Wei [4]. They considered homogeneous 
and inhomogeneous processes, and a Choleski decomposition of the covariance 
matrix was used to correlate the Monte Carlo sampling. The direct sampling 
method, coupled with the finite element procedure were used to study porous 
random fields for groundwater flow. 

The Neumann expansion technique proposed by Shinozuka et al. [5,6] 
used sampling techniques that successfully reduced the computational effort in 
solving the finite element equations independently for each sample. The 
Neumann expansion method effectively employs a perturbation expansion in 
conjunction with the Monte Carlo simulation so that only a single stiffness 
matrix factorization is required. The method allows for large variability of the 
random variables to be modeled without loss of accuracy. It can be easily 
adapted to an existing finite element code with little change. 

Contreras [7] proposed a different method in which stochastic differential 
and difference theory is applied to structures discretized using the finite element 
method. A semi— discretized formulation is employed in which a finite 
difference method is used for the time domain and the finite element method is 
used for the spatial domain. This technique involves a complete reformulation 
of an existing finite element program in order to model the stochastic 
differential equations. 
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A method known as stochastic linearization [8,9] has been used to solve 
nonlinear stochastic differential equations of random vibration in which the 
loading is the only random process. Recently, Mohammadi and Amin [10] and 
Casciati and Faravelli [11—15] employed it for dynamic analysis of material 
nonlinear continua. The stochastic linearization technique is used to linearize 
the nonlinear hysteresis behavior in the model. In addition, Chen and Yang [16] 
used a finite element formulation combined with stochastic linearization and 
normal mode methods to study geometrically nonlinear random vibration of 
plate and shell structures. In general, this method has only been applied to 
random loading problems, and not to problems in which uncertainties exist in 
the properties. 

Faravelli [17] has demonstrated another statistical approach in which a 
planned set of experiments around the space of the central values of the 
different random vectors is performed using a standard finite element code to 
determine the response to the different inputs. A regression analysis is then 
used to fit the response to an appropriate polynomial of the input variables. 
A level— 2 reliability approach is introduced to obtain approximations of the 
cumulative distribution functions of the response variables. Nonlinear problems 
can be solved as well as linear. In order to minimize the number of numerical 
experiments necessary, experimental design theory is used. It seems that this 
method would be computationally costly when a large number of random 
variables is involved. In reference [18], the method was applied to an 
automotive impact problem. A similar method was used by Chryssanthopoulos 
et al. [19] to perform a reliability— based design of stringer stiffened cylinders 
under axial compression. 
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Larder [20] developed a method for studying failure in parallel fiber 
composites. A deterministic finite element formulation was used to model the 
fiber and matrix of a representative cell, and then Monte Carlo simulation 
methods were used to simulate randomly occurring damage. The program 
however, could not model general structural level problems. 

Deodatis and Shinozuka [21] developed a probabilistic model for the 
spatial strength variation in laminated orthotropic composites. Monte Carlo 
simulation techniques were coupled with Tsai— Hill and Tsai— Wu failure criteria 
and assumed failure mechanisms. 

Finally, Ditlevsen et al. [22] have conducted research in improving the 
efficiency of the Monte Carlo method. An off-mean centered directional 
importance sampling procedure is compared to other uniform directional 
sampling methods. The aim was to further reduce the number of samples 
required as compared to other improved Monte Carlo procedures. 

2.1.2 Non— Statistical Approaches 

Non-statistical finite element approaches seem to be getting the most 
attention in current research. These approaches include second— moment 
analysis, numerical integration, and the new iterative perturbation methods, 
which use multiple regression as a post— processor. These techniques will be 
discussed and compared in detail in the following review. 

Ang and Tang [23], discussed a method of using the Taylor series 
expansion of a general function g(X) about the mean of the random variable X 
and truncating the series after the second order terms. They showed that the 
mean of g(X) can be estimated using a second-order approximation which uses 
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the first and second statistical moments of X, and that the variance of g(X) can 
be estimated using a first-order approximation which requires only the second 
moment of X. Since higher statistical moments of the original variate X would 
be required for a higher-order approximation to the variance of g(X), the 
variance was left at first— order. Higher moments of the original variates are 
generally not known. A major advantage of this method is that the 
multivariate distribution function does not need to be known, but only the first 
two moments. Since a first-order Taylor series approximation is used for the 
variance, then uncertainties in the original variates cannot be too large. This 
means that deviations from the mean of the random variables of the function 
cannot be large, typically not greater than 10 percent. However, they showed 
that acceptable results can still be obtained for coefficients of variation as high 
as 20 percent. An important point is that these approximations have proven to 
be adequate even when the function g(X) is nonlinear, as long as the variance of 
X is small relative to g(X). This leads directly to the works by Liu et al. 
[33-40] to be discussed later. 

Handa and Andersson [24] used the above ideas combined with the finite 
element method to study linear truss structures and the effects of correlation. 
At about the same time Hisada and Nakagiri [25] used the same "perturbation" 
(second— moment) method to analyze structures with uncertain shapes. In this 
work the geometrical coordinates were allowed to be random variables. 
Nakagiri et al. [26—29] were the first, and to the author’s knowledge, the only 
researchers to apply this method to the analysis of composites. In [26], 
probabilistic eigenvalue analysis of linear vibration of fiber reinforced laminated 
plates was studied, in which the fiber orientation angle and the layer thickness 
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were random variables. The effects of these random variables on the 
eigenvalues of the plate with and without correlation was investigated. In [27], 
reliability indices of a fiber reinforced laminated plate were calculated as the 
end result of the stochastic finite element technique. Here elastic constants of 
the layers were treated as the random variables, but no spatial correlation was 
considered. The strengths used in the failure criterion to determine reliability, 
were also variates. Cases in which flat plates under uniaxial tension, with and 
without holes, were analyzed to investigate the effect of the distributed stress 
on the reliability index. In [28], Tani and Nakagiri extended the method to 
perform design optimization of the reliability index of fiber reinforced plastic 
(FRP) laminated plates with probabilistic elastic constants and lamina 
strengths. The fiber orientation angle was the design parameter in the 
optimization, and failure of each lamina was determined by the Tsai— Hill 
criterion. Sato, Watanabe, and Nakagiri [29] used the second— moment method 
to perform a reliability analysis of a FRP pressure vessel with probabilistically 
distributed stacking parameters. Linear elastic material and geometrically 
linear behavior was studied in their work. 

A different approach to the second-moment method was used by 
Lawrence [30]. Here expansions of random variables in terms of orthogonal 
Galerkin type functions became the core of his "basis random variable" 
technique. A Rayleigh— Ritz finite element model was derived, and the first and 
second moment statistics were obtained from a series of linear solutions. 
Lawrence [30] achieved the same accuracy as other second-moment 
FEM techniques; however, he stated that the greatest advantage of his method 
over others is the ease of application. Using the basis random variable 
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representation, the random character of the problem is treated "merely as an 
extra set of dimensions". Later, Lawrence [31] applied the method to 
probability based design, in which sensitivities of reliability indices to design 
variables are calculated to enhance the design process. Additionally, Ghanem 
and Spanos [32] developed an independent but similar Galerkin— based response 
surface approach. In these works only linear problems were studied. 

Liu, Belytschko, and Mani [33-40] have carried out the most 
developmental work using the second— moment finite element method. They 
have applied the method to nonlinear structural dynamics for correlated or 
uncorrelated discrete random variables [34,35], and later to nonlinear structural 
dynamics problems with both homogeneous or inhomogeneous random fields 
[33]. They termed the method the Probabilistic Finite Element Method 
(PFEM), and applied it to geometric as well as material nonlinear problems. 
Good results were achieved since the coefficients of variations were kept around 
10 percent. It is noted that the probabilistic density functions must also be 
considered to have decaying tails to maintain good results. Only isotropic 
materials were studied. Two methods for improving computational efficiency 
were developed: 1) transforming the correlated random variables to independent 
uncorrelated random variables thereby reducing the order of the tensor 
multiplications required, and 2) an adjoint method was introduced to reduce 
computations by allowing the statistics of the response to be calculated only in 
a reduced domain from the original problem. 

Liu, Besterfield, and Belytschko [36] developed a probabilistic 
Hu— Washizu variational principle (PHWVP) for the PFEM. They applied it to 
both linear and nonlinear elasticity. Using this principle, the probabilistic 
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distributions for the constitutive law, compatibility condition, equilibrium, 
domain and boundary conditions can be incorporated into the PFEM. Thus all 
aspects of the problem were treated as random variables or fields. 

The latest work by these authors involves an optimum fusion of 
PFEM and reliability analysis. Liu, et al. [41] applied the method to 
probabilistic fracture mechanics. Uncertainties in loads, material parameters, 
geometry, crack length, orientation and location can be incorporated into the 
problem to determine the probability of brittle fracture. The reliability analysis 
is performed using a Kuhn— Tucker optimization procedure. An enriched finite 
element with the embedded crack— tip singularity is employed in this overall 
reliability package. In [42], they extended the method to fatigue crack growth 
problems. Liu, Chen, and Lu [43] have also applied the reliability and 
PFEM fusion method to the fluid— structure interaction problem in which the 
influences of random parameters on the energy transfer between the structural 
system and the acoustic system were studied. 

Several other authors have fused the first and second— order reliability 
methods with finite element methods. Here the differentiation of the finite 
element equations with respect to the random variables is performed as in the 
second— moment formulation, except that these are then blended directly into 
the first-order reliability equations. In this way the small probabilities of 
failure required in structural reliability calculations can be accurately predicted. 
These small probabilities at the tails of the distributions are poorly defined by 
first and second moments. Ambjerg— Nielson and Bjerager [44] applied a form 
of the method to linear and nonlinear truss structures. Der Kiureghian et al. 
[45,46] applied it to linear and nonlinear beam and plate problems. The key 
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point is that while the PFEM methods evaluated statistical moments, this step 
is avoided altogether here in the process of estimating reliability. 

Other non— statistical techniques have been developed and are discussed 
next. Liu et al. [34] developed a numerical integration procedure in which 
Hermite— Gauss quadrature is used for numerically evaluating the integrals 
involved in the definition of expectation and variance. The method assumes 
uncorrelated fields and has limited use, especially for problems involving large 
numbers of random variables (also see Gorman [47]). Takada and Shinozuka 
[48] have developed a new method in which local integrations of the continuous 
stochastic field are made on an element by element basis to form element 
stiffness matrices so that the random field is transformed into only a few 
random variables. This new method has been proven to allow much coarser 
meshes than the first-order second-moment methods, but requires more 
computational time for equivalent mesh sizes. It has only been applied to linear 
elastic problems. A similar idea is used by Weiqiu and Weiqiang [49], in which 
local averages of the random field are used to improve the accuracy of the 
second-order perturbation method. They applied it to random eigenvalue 
problems. 

The best alternative to the first-order second— moment methods is one 
developed by Dias and Nagtegaal [50]. They introduced stochastic finite 
element techniques based on the Fast Probability Integration (FPI) algorithms 
developed by Wu and others [81—85]. Using FPI, and given a closed— form 
expression of the response variable as a function of the uncorrelated random 
variables (normalized to zero mean and unit standard deviation), an 
optimization algorithm is then used to determine the minimum distance from 
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the origin to a given limit state. This minimum distance is commonly referred 
to as the safety index, from which the probability of exceedance of the limit 
state can easily be computed. Thus what is needed is a method to relate the 
response of the structure to the uncorrelated random variables. By performing 
a spectral decomposition of the covariance matrix, the correlated random 
variables can be related to a set of uncorrelated random variables. In order to 
determine the response of the structure, the deterministic finite element model 
is used in conjunction with a prescribed set of perturbations of the input 
random variables. To obtain a closed— form expression of the response, a 
multiple regression procedure that fits the stored perturbations to a polynomial 
expression is used. 

Obviously the method requires a multiple set of solutions of the finite 
element model which can become computationally expensive. In order to 
overcome this problem, an approach based on a modified Newton iteration 
method was developed. Basically it uses an iterative perturbation procedure 
about the known mean deterministic solution, so that multiple perturbation 
solutions can be obtained using only one factorization of the stiffness matrix. A 
notable advantage of this method is the lack of modifications required to 
incorporate it into an existing finite element program. In addition, the method 
is not limited to any specific probability distributions for the input random 
variables. What is needed is an efficient data storage structure to keep track of 
the multitude of perturbation solutions. Many other authors have worked with 
Dias in this arena (see [50—55]). It is of importance to mention that the 
Probabilistic Structural Analysis Methods (PSAM) being developed at 
NASA— Lewis Research Center use this technique, due to its ease of 
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incorporation into any finite element formulation by simply coupling it with 
existing modules for reliability computations. 

Various examples exist in the literature of application of the previously 
discussed methods to specific problems. A recent book edited by Liu and 
Belytschko [56], contains a wealth of interesting examples. Recent works using 
the first-order second-moment stochastic FEM method for analysis of soil 
structures are given in references [57,58]. Various stochastic FEM methods for 
solving material nonlinear problems are reviewed in references [59—62] . 

2.1.3 Comparison of Approaches 

From the literature review, it was quickly determined that the 
non— statistical techniques would prove to be the most promising for application 
to linear and nonlinear composite structures. It was also noticed immediately 
that little work has been done (only Nakagiri [26—29]) in the area of stochastic 
finite element procedures for composite structures. Out of all the works 
reviewed, the second-moment techniques have received the most attention for 
nonlinear applications, but only for isotropic behavior (see Liu et al. [33^0] 
and Liu and Der Kiureghian [46]). The iterative perturbation techniques also 
seem very promising. The following discussion is an attempt to compare the 
pros and cons of these two procedures. 

The first-order second-moment approach is basically a mean-centered 
perturbation method since the Taylor series expansion is performed about the 
mean point. Since it is only a first-order approximation for the covariance, 
then the method does not allow for large variations of the input random 
variables from the mean. Typically coefficients of variation of 10 percent give 
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good results, however, and for most material property variations this range is 
well within the expected level of uncertainty. 

In order to apply the second— moment technique to an existing finite 
element code, significant changes must be made in order to calculate partial 
derivatives of the stiffness matrices and loads with respect to the random 
variables. These derivatives can be obtained either using analytical 
differentiation of the formulation under consideration or finite difference 
techniques. The latter method should allow more generic programming 
techniques to be employed, making it more formulation independent. As the 
number of random variables in the problems increase, it becomes apparent that 
a significant amount of storage becomes necessary to retain all the partial 
derivatives. 

An advantage of the second— moment method is that it has been 
previously demonstrated for nonlinear problems by Liu et al. [33 — 40], with good 
success. The authors also developed techniques for efficiency improvements to 
reduce the computational costs considerably. 

Once the statistical moments of the response are determined using the 
second— moment formulation, any of the reliability estimation techniques can be 
employed. Recent work however, has shown that this is not necessary; in fact, 
it is more accurate not to calculate the moments but to evaluate the sensitivity 
derivatives in the finite element program and use them directly to estimate 
reliability. In any case, unlike the iterative perturbation methods, no regression 
curve fitting is required in this step. 

The iterative perturbation methods involve a higher order perturbation 
than the strict first— order second— moment methods. Using the iterative 
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procedure, higher order terms are effectively captured. The stability of the 
algorithm is conditional however, and convergence is not assured for large 
stiffness perturbations. The method can be generalized so that any finite 
element formulation can be incorporated into the algorithm without any 
element level program changes. Like the previous second— moment methods, it 
requires significant storage in order to retain all the response solutions to the 
various perturbations of the random variables. This storage coordination must 
be very well organized and self contained as well. 

The iterative methods have just recently been demonstrated for material 
nonlinear problems in the current literature. Convergence stability problems 
have been reported when constraint equations exist such as those for thin plates 
and shells allowing transverse shear deformation, deviatoric rate— independent 
plasticity, and incompressible or near— incompressible elastic materials. 
Currently, stability limits are being developed to alleviate these problems. 

The results of the iterative perturbation finite element program are data 
files of response for various perturbed inputs. In order to make reliability 
calculations, the FPI reliability algorithms require a closed form equation 
involving the response variables as a function of the input random variables. 
Thus multiple regression techniques are employed at this step to obtain a 
polynomial equation from the response data files. 

After reviewing all of the above facets of the two methods, it was decided 
to use the second— moment methods for the present application. Since an 
experience base exists for the nonlinear applications, and since large 
non— normal variations are not anticipated for the uncertainties, the methods 
should be acceptably accurate. As for the storage problems, it appears that 
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both methods will require considerable storage, so this really is not a factor. 
Since Liu et al. [33-40] has developed techniques for reducing computational 
costs, it was anticipated that the effort spent making significant programming 
changes would be rewarded by a reduction in CPU charges. This is very 
important since a rather large number of random variables are needed to model 
the uncertainties in a composite laminate. 

2.2 Micromechanics Models 

One of the goals of this research is to model from a micromechanics level 
the uncertainties in laminated composites, involving both linear and nonlinear 
behavior. It was discovered, however that certain algorithmic difficulties 
existed in coupling a micromechanics— based plasticity theory with the 
second— moment stochastic finite element methods. Thus the micromechanics 
model was limited to linear and geometric nonlinear elastic behavior for this 
study with the possibility of adding the material nonlinear case at a later date. 
In any case, it was desired to select a micromechanics theory that has proven to 
be effective in the literature. Arenburg [63] has given a very thorough review of 
this area of research, and has recommended the model by Aboudi [64], This 
model is based on a higher-order continuum theory that can account for 
particulate or continuous fiber reinforcement, general mechanical and thermal 
load histories, and damage in the form of fiber— matrix debonding. Arenburg 
used this model to develop a finite element program for the analysis of 
metal— matrix composite plates. In his work, the unified theory of Bodner [65] 
was used to model the nonlinear viscoplastic behavior of the matrix. In the 
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present work, Aboudi’s model will again be used, but in this case will be limited 
to elastic behavior. 

2.3 Macroscopic Anisotropic Plasticity 

In order to model the nonlinear material behavior existing in metal 
matrix composites, it was decided to use the macroscopic or smeared approach. 
This method proved to be more compatible with the probabilistic finite element 
method both from an algorithmic and a computational expense point of view. 
Many authors have made contributions in this area. Griffin [66] and Arenburg 
[63] have presented surveys of the theories of plasticity in anisotropic materials. 
It is not of interest here to review all of these, but to touch on the most notable 
and then explain the reasons for selecting the theory used in the present work. 

Hill’s anisotropic theory of plasticity [67,68] has received much 
attention. This theory was based on a generalization of the von Mises yield 
criterion which assumed yielding was independent of hydrostatic stress and that 
plastic flow was incompressible. Whereas these assumptions are standard and 
experimentally validated for metals, they have been shown to be incorrect for 
some materials by Lin et al. [69]. Griffin [66] and Chandrashekhara [70] have 
applied the theory in finite element formulations. While Griffin showed 
generally good agreement with off-axis test data, Chandrashekhara did not 
make any experimental comparisons. 

An endochronic theory of plasticity of Valanis [71] was extended to 
transversely isotropic media by Pindera and Herakovich [72]. Various effects 
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such as cyclic loading and unloading and stress interaction were successfully 
demonstrated. 

Recently, Sun et al. [74—77], presented a plasticity formulation based on 
the anisotropic plastic flow rule proposed by Kachanov [73]. The yield function 
they selected was quadratic in stresses and, in general, excludes the 
assumptions of incompressibility of plastic strains and that no yielding is caused 
by hydrostatic stresses. The formulation was installed in a finite element 
program, and much experimental validation work was performed for 
Boron/ Aluminum composites [74—77]. Examples of off-axis uniaxial test 
specimens with holes, edge cracked panels, and cyclic loading and unloading 
under constrained plasticity conditions were demonstrated with good 
analytical/experimental comparisons. Due to the generality and amount of 
experimental work done by Sun and his coworkers, this plasticity formulation 
was selected for the present work. Algorithmic modifications were made for 
computational efficiency and compatibility with the second— moment 
probabilistic method. 

2.4 Reliability Estimation 

Orginally it was planned to estimate reliability as a post— processing 
function using the first and second moments of the stress or displacement 
response output from the probabilistic finite element analysis. However recent 
work [41-42,45-46], which fused the second— moment probabilistic finite 
element methods with the first and second— order reliability methods, have 
shown that direct use of the response derivatives with respect to the random 
variables (obtained from the finite element program) in the reliability 
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algorithms is much more accurate for small probabilities of failure. These small 
probabilities at the tails of the distributions are poorly defined by the first and 
second moments. For this reason, the reliability estimations in the present 
work are limited to approximations based on the first and second moments 
(these moments are still quite useful for sensitivity analysis). Future work will 
use the response derivatives already evaluated in the finite element program 
directly for reliability estimation. Appendix A gives a brief explanation of the 
reliability procedures. 

In reviewing the vast literature in reliability computation, the first and 
second— order reliability methods are prominant. References [23,45,46,56] give a 
good description of these methods. The family of Fast Probability Integration 
(FPI) algorithms has also been found to be important in recent work. The 
origin of the FPI methods can be traced to the definition of a safety index 
proposed by Hasofer and Lind [78]. Rackwitz and Fiessler [79] later extended 
the algorithm to accommodate non— Gaussian distributions. Recent work by 
Wu et al. [81—85] contains further improvements in terms of accuracy and 
computational efficiency. 
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3. FINITE ELEMENT FORMULATION 


In this chapter a review of the finite element formulation of the 
incremental equations of motion of a continuous medium including geometric 
nonlinearity will be presented. The degenerated 3— D shell element for 
composite laminates will then be reviewed. The formulation is based on the 
work by Liao and Reddy [86]. Summation convention on repeated indices is 
assumed throughout this study. 

3.1 Principle of Virtual Displacements 

The principle of virtual displacements requires that 



(3.1) 


where 6 is the variational symbol 



the Cartesian components of the Cauchy stress tensor at 
time t + At (i.e., configuration 2). 



the Cartesian components of the infinitesimal strain 
associated with the incremental displacements Au- in 
going from the configuration at time t to the configuration 


at time t + At, i.e. e j; = j ( y - + a x which are also 

J j i 

referred to this unknown configuration at time t + At. 
Cartesian components of a point in the configuration at 
time t + At, 


<9Am 5Au- 


and 



(3.2) 
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In equation (3.2), t^ and are the components of the externally applied surface 
and body force vectors respectively, and is a virtual variation in the current 
displacement components u^, where 


being the Cartesian coordinates of a point in the configuration at time t = 0. 

3.2 Total Lagrangian Formulation 

Equation (3.1) cannot be solved directly since the configuration at time 
t + At is unknown. An approximate solution of equation (3.1) can be obtained 
by referring all variables to a previously calculated known equilibrium 
configuration and linearizing the resulting equation. The total Lagrangian 
(T.L.) formulation, in which all static and kinematic variables are referred to 
the initial configuration at time 0 of the body, is used here. The applied forces 
in equation (3.1) are evaluated using 


where T^, F^, dA, dV, p Q refer to time t = 0 and t^, f^, da, dv, p refer to time 


The volume integral of Cauchy stresses times variations in infinitesimal 
strains in equation (3.1) is transformed to an integral over a known volume (the 
initial volume) using the 2nd Piola— Kirchhoff stress tensor and the energetically 
conjugate Green— Lagrange strain tensor to give 



T^dA = t^da 

W v = 'V 1 ' 


(3.3) 


t -f- At. 




where 
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Cartesian components of the 2nd Piola— Kirchhoff stress 
tensor corresponding to configuration at time t + At but 
measured in the configuration at time t = 0. 

Cartesian components of the Green— Lagrange strain tensor 
in the configuration at time t + At, referred to the 
configuration at time 0: 


du- 


E • = i (u- • + u. . + u, u, .) and u- . = 
ij 2 ' i,j j,i k,i k,j' i,j 3XT 


The 2nd Piola— Kirchhoff stress tensor referred to the configuration at time 
t = 0 is defined as 


S 


ij 


rVsr X j,r 


(3.5) 


Substituting equations (3.3) and (3.4) into (3.1) and (3.2), we obtain 


J V W V = “ 


(3.6) 


which involves the equilibrium for the body in the configuration at time t + At 
but referred to the configuration at time t = 0. In equation (3.6), R is 
calculated using 


R = 


T kV A 


+ 


W*k dV 


(3.7) 
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The stresses S-- and strains E-., 
^ J j 


decomposed as 


which are unknown at time t + At, are 


S.. = S 1 . + S°. 
U iJ iJ 


E.. = E*. + E°. 
iJ iJ 


(3.8) 

(3.9) 


where S*. and E*. are the known 2nd Piola— Kirchhoff stresses and 
^ J ^ J 

Green— Lagrange strains in the configuration at time t, and S?j and E?j are the 
incremental components of the same quantities at time t + At. It follows from 
equation (3.9) that 


E?. = e? - + v°- 
iJ U 'iJ 


where 


(3.10) 


e °j - 5 (Au i,j + Au j,i + U k,i Au k,j + Au k,i U k,j> 


(3.11) 


= linear part of strain increment E?. (linear in u- •) 

* J 1 > J 


"“j = \ ( Au k,i Au k,j) 


(3-12) 


= nonlinear part of strain increment E°j 


and 


26 


dU ■ <9Au, 

U i,j = 3XT ’ Au k,i = ~dT~ 

.j j i 

In the above formulation, the definition of the Green— Lagrange strain tensor 
was employed along with the following definition 

u i = u i + Au i 

where Au- is the incremental displacement and U- is the displacement at time t. 

The constitutive relationship for the incremental 2nd Piola— Kirchhoff 
stresses and incremental Green— Lagrange strains is assumed to be governed by 
the generalized Hooke’s law, 



= C. • E 
ljrs rs 


(3,13) 


where C-. are components of the constitutive tensor. Substituting equations 
ljrs 

(3.8)— (3.13) into equation (3.6) it follows that 


L w>V v + L s u<j dv = R - L s ij fe ?j dV < 314 ) 


which represents a nonlinear equilibrium equation for the incremental 
displacements Au^. 


3.3 Linearization of Incremental Equations of Motion 

Since equation (3.14) is nonlinear in the incremental displacements and 
cannot be solved directly, approximate solutions are obtained by assuming 
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E?: = e? : and <5E? : 


ij iJ 


ij 


<5e j j. The constitutive equation becomes 


S°. = C-. e° 
ij ijrs rs 


(3.15) 


Making the above substitutions the approximate equation is now 


V 


C-. e°>?.dV + f S; .fy°.dV = R - [ s] .fe°.dV 
ijrs rs ij J v ij ij J v *J *J 


(3.16) 


Using Hamilton’s principle, the effects of the inertial forces can be included. 
Employing similar procedures as before, we obtain the equations of motion of 
the moving body at time t + At in the variational form as 


, vW v + j v + J v s ^ij dv = R - | v s u*V v 

(3.17) 

Equation (3.17) forms the theoretical basis for the finite element model. 


3.4 Finite Element Model 

Using standard isoparametric interpolation, the final incremental 
equations of motion for an element are given by 

[M e ]{A e } + ([K«l + [K« L ]){A e } = {R e } - {F e } (3.18) 

p 

where {A } is the vector of nodal incremental displacements from time t to 
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time t + At in an element, and [M]{A}, [KjJ{A} ( [K^^]{A}, and {F} are 
obtained, respectively, from the following integrals 


V 


/> 0 V u i dV > 


C-. e° <5e°.dV 
y ijrs rs ij 


V 


sfj^jdV , and 


Q-*- fie 0 dV 
v ij i J V 


These integrals can be written in the standard matrix form as [86] 


[K L ] = . 

v [B L ] T (C][B L ]dV 

(3.19) 

S! 

tr 1 

II 

v ^ B NL^ S ^ B NL^ dV 

(3.20) 

[M] = 

[ P 0 !H] T [H]dV 

(3.21) 

{F} = 

|B L ] T {S}dV 
Jy L 

(3.22) 


In equations (3.19)— (3.22), [BjJ and [B^jJ are linear and non-linear 
strain-displacement transformation matrices, [C] is the incremental 
stress— strain material property matrix, [S] is a matrix of 2nd Piola— Kirchhoff 
stress components, {S} is a vector of these stresses and [H] is the incremental 
displacement interpolation matrix. All matrix elements correspond to the 
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configuration at time t and are defined with respect to the configuration at 
time 0. 

The finite element equations in (3.18) are second-order differential 
equations in time. Equation (3.18) must be further discretized in time to obtain 
algebraic equations, which can be assembled and solved after imposing initial 
and boundary conditions. Here the Newmark integration scheme is used to 
approximate the time derivatives. The resulting algebraic equations are given 
by (see Reddy [87]) 

[K]{&} = {R} (3.23) 

where {A} is the vector of nodal incremental displacements at time t, and 


[K] = aJM] + [K l ] + [K nl ] 


{R} = {R}-{F} + [M](a 1 {A} + a 2 {A}) 


a. ^ , a, - a Q At , a 2 - 1 


«At)‘ 


(3.24) 


where a = 1/2, ft = 1/4 for the constant average acceleration method and At is 
the time step. Once equation (3.23) is solved for {A} at time t + At, the 
acceleration and velocity vectors can be computed from 


{ 2 A} = a 0 { 2 a}-a 1 { 1 A} -a/A} 

{ 2 A} = {*A> + a 3 { 1 A> + a 4 { 2 A} (3.25) 
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where the superscript 2 refers to time t + At and 1 refers to time t, and 

a^ = (1 - a)At , a^ = aAt. 

Equations (3.23) and (3.24) are only a linearized version of the actual 
governing equations of motion. Hence equation (3.23) should be solved 
iteratively for each time step until the actual equations of motion are satisfied 
to a required tolerance. Here the Newton— Raphson iteration technique is 
employed. In the Newton— Raphson method, the equation to be solved at the 
ith iteration for time step t + At has the form 

(a 0 (M] + [K l ] + [K NL ])(‘ -1 )({A}W - {A}* 1 - 1 )) 

= {R} - (M](a 0 {A>( i - 1 > - aj{ A} - a 2 { A}) - {F}^ i— ^ (3.26) 

where 

{F}( 1—1 ) = J [B^ 1) ] T {S}^ 1—1 ^dV 

After the ith iteration, the nodal displacement at time t + At is updated by 
{ 2 A }(i) = { 2 A}( i_1 ) + ({A}^ - {A}( i_1 )) 

in which 

{A} (0) = {0} 

3.5 Degenerated Laminated 3— D Shell Element 

The degenerated shell element developed by Chao and Reddy [88] for 
composite laminates will be utilized here. The shell element is degenerated 
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from the three-dimensional solid element by imposing two constraints: (1) 

straight lines normal to the midsurface before deformation remain straight but 
not normal after deformation; (2) the transverse normal components of strain 
and hence stress are ignored in the development. The resulting non-linear 
formulation allows arbitrarily large displacements and rotations of the shell 
element but small strains, since the thickness does not change and the normal 
does not distort. 

As in [86,88], the 3— D solid element in Figure 3.1 is the starting point. 
The curvilinear coordinates in the middle surface of the shell are given by (, tj, 
£, which are normalized such that they vary between —1 and 1. Here £ denotes 
the thickness coordinate. The coordinates of a point are given by 



n 

S 

k=l 



(3.27) 


where n denotes the number of nodes per element, and <p, (£,rj) are the 

Lagrangian finite element interpolation functions. The normal vectors 
connecting the upper and lower surfaces of the shell at node k are defined by 


their components as 


V 3, = ( X Xp- (^bottom 

(3.28a) 

4 = ^/|V^| 

Substituting these expressions into (3.27), we obtain 

(3.28b) 
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Figure 3.1 General 3-D solid geometry and resultant shell element. 
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(3.29) 


X i = + Kil 

where h^ = |Vg| is the thickness of the shell at node k. The displacement 
components become 

U i = % - X i = k l 1 W-M + i h k(‘ e 3i - °«Si» ( 3 3 °) 

Auj = - U, = * k «,»)[Au* + $ h k ( t+A ^ - *4,)] 

(3.31) 

t t k o k 

where x- are the Cartesian coordinates at time t and eg- and e 3i are the unit 

vectors at time t and 0, respectively. In order to update the vectors in (3.28) 
for small rotations, the rotation dfl is expressed as 

dS = ^ ^ 

t k 

The increment of eg can be written as 

= dfi * ‘ejj = ^ t e \ - ^ ^ (3.32) 

The equation (3.31) becomes finally 
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Au i= £ + (3.33a) 


k=l 


t k t k 

where e.^ and eg are determined from 


1 ip t k. 

| Eg x e 3 1 

^ = ^3 * l e\ (3.33b) 


E- are the unit vectors of the stationary global coordinate system. Equation 
(3.33a) can be expressed in matrix form as 


{Au} = {Auj^ Aug = [H]{A e } 

where 

{A e } = {Au^ Aug Au^ 0^ 0^}^ (k = 1 to n) 


and [H] is the incremental displacement interpolation matrix, which can be 
found in [ 86 ]. 

Next, the linear strain increments {e 0 } = {e°^ e^g e^ 2e°g 2 e °3 Se^}^ 
are written in matrix form 


{e 0 } = [A]{ Au} 


where 


{Au} = {Au,, Au, 2 Au 13 Au 2i , Au 2 2 Au^ Au^, Au 32 Auj j} 
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and 


U i>j 


auj 

3X- 


The matrix [A] can be found in [86]. Denoting {Au} by 

{Au} = [N]{Au} = lN][H]{A e } 

with [N] being a matrix of differential operators, then {e 0 } becomes 

{e 0 } = [A]{Au} = [A][N](H]{ A e } i [B L ){A e } 

The usual Jacobian transformations are used to express u- • in terms of the 

*>J 

derivatives with respect to the global coordinates. Also in this work the 

/ / 

integration is done in terms of the local curvilinear coordinate system (xp x 2 , 

/ 

Xg), as appropriate transformations are made. 

/ 

The constitutive matrix [C ] for a lamina in the local coordinate system 
is given by 



where 


c n 

C 12 

C 13 

0 

0 


C 12 

C 22 

C 23 

0 

0 


C 13 

C 23 

C 33 

0 

0 


0 

0 

0 

C 44 

C 4S 


0 

0 

0 

C 45 

C 55_ 


C ll = 

= m 4 Q 

n + 2m 2 n 2 (Q 

12 + 2( ^33^ + n ^22 

C 12 = 

= mV 

! (Q u 

+ ^22 

4 «33 

) + (m 4 + n 4 )Q 12 


(3.34) 
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C 13 = mn[m 2 Q n - n 2 Q 22 - (m 2 - n 2 )(Q 12 + 2Q 33 )] 
C 22 = n 4 Q n + 2m 2 n 2 (Q 12 + 2Q 33 ) + m 4 Q 22 
C 23 = mn[n 2 Q u - m 2 Q 22 + (m 2 - n 2 )(Q 12 + 2Q 33 )] 
C33 = m n (Q n + Q 22 - 2Q 12 ) + (m - n ) Q 33 

C 44 = ^44 + n ^55 ’ C 45 = mn ^44 - Q55) 

/ 2 2 

C 55 = m Q55 + n Q44 » m = cos 6 , n = sin 6 


and Q-j are the plane stress— reduced elastic coefficients in the material 
coordinates. By neglecting the normal stress in the thickness direction, the 
stiffnesses Q-j for an orthotropic lamina are determined from the 
three-dimensional orthotropic stiffnesses, which can be expressed in terms of 
engineering constants as, 


E 


Qll — 4 _ jy ^ > Q 

11 1 ^ 12^21 


12 


^ 12 E 22 
1 - ^ 12^21 


^22 “ I 


E 22 

v \2 u 2l 

( 3 . 35 ) 
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Q 44 ~ ( G 13^ K ’ ^55 “ ( G 23^ K ’ ^33 “ G 12 


where K is the shear correction coefficient taken to be 5/6. 

In general, Gauss quadrature is used to integrate the element matrices. 
For laminated composite structures, the thickness direction integration could be 
done using Gauss quadrature, or the problem could be reduced to a 2— D one 
and hence explicit integration performed. In order to do this, the Jacobian 
matrix dependency on the thickness coordinate ( is neglected. This is a valid 
assumption if the thickness to radius of the shell ratios are small. In this study 
we use the explicit integration for elastic behavior and numerical integration 
(quadrature) when elastic— plastic behavior is assumed. 

From [86], the following symbolic representations of desired quantities 
are expressed explicitly in terms of 

{Afl'} = ([DHl] + <[DH2]){A e } 

[A] = [SD] + C[TD] 

[B l ] = ([SD] + C[TD])([DH1] + ([DH2]) 

(3.36) 

= [SD1] + C[SD2] + < 2 [SD3] 

{Ejj} = {S^ + <{S2} + C 2 {S3} 
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[B nl ] = [DH1] + <[DH2] 

These representations can be used to perform explicit integration through the 
thickness. 

In order to avoid locking behavior for thin shell structures, selective 
reduced or fully reduced integration is employed. Deformation dependent 
loading is also used and is described in [86]. 
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4. MICROMECHANICS 


4.1 Introduction 

The Aboudi theory has been included in order to study the effect of using 
micromechanics constituent properties as random variables. In this section a 
review of Aboudi’s theory is discussed, which is modeled after the development 
in [64]. 


4.2 The Aboudi Micromechanics Model 

This theory involves the solution of a suitable boundary value problem 
whose domain is a typical representative volume V. The composite is modeled 
as an isotropic viscoplastic matrix reinforced by an elastic transversely isotropic 
fiber of rectangular cross section. The fibers extend in the x^ direction and are 
arranged in a doubly periodic array in the X 2 and x^ directions as shown in 
Figure 4.1a. The rectangular fiber has cross sectional dimensions hp with 
1^, ^ denoting the matrix spacing. Figure 4.1b shows the representative cell 
necessary for analysis due to the periodic arrangement. The cell is further 
divided into four subcells a, 0 = 1,2 each with a local coordinate system (xp 

*(«) xM) 

*2 ’3 '' 

The displacement field considered by Aboudi is a first-order expansion 
in each subcell, since only average behavior of the composite is sought. This 
displacement field of Aboudi is given by 


= + + , i = i )2 ,3 (4.1) 


where wj°fl are the displacement components of the center of the subcell and 
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3 


(b) The representative cell. 


Figure 4.1 Micromechanics subcell geometry 
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<py*P' and V'j ^ characterize the linear variations of the displacements within 
the subcell in the and directions, respectively. In this section, 

repeated a and 0 do not imply summation. 

Displacement continuity at the interfaces between the subcells are given 

by 


w( U > = w! 12 > = w( 21 > = w! 22 ) , W. 
11111 

(4.2) 

h l 4 1 ® + = ( h l + h 2^ od; 

(4.3) 

l \ 4 al) + hA a2) = ^1 + V ^ 

(4.4) 


The average strain in the representative cell is given by 


£ 


ij 


1 


2 

E 


V a, 0=1 



(4.5) 


where v^ = h^ and V = (1^ + h 2 )(^ + 4,) is the area of the representative 
cell. The infinitesimal strain tensor is 


1 J z 1 J J 1 


ij = 1,2,3 


(4.6) 


where 


d 


d 


d 


Combining equations (4.3), (4.4), (4.6) then the following result is obtained 


, <9w. dw- 

~ e ii = 2 + 33$ 

J J i 


(4.7) 


The average stress in the composite cell is 


a- ■ — — £ v rSi° 

V V#1 




(4.8) 


where s| ^ is the average stress in the subcell. sj is determined by 


( 03 ) _ l r h “/ 2 [ l M 2 

ij “ V J J ‘J 2 3 


(4.9) 


h a /2 ^/?/2 


Aboudi’s theory is designed to allow both fiber and matrix constituents 
to be elasto— plastic materials. Thus the strain rate in the subcells are 
decomposed into elastic and plastic parts as 


.(a/3) = -E(a/3) -P(a/3) 

ij ij ij 


(4.10) 


The constitutive equation for transversely isotropic constituents in the 
elastic range, with x^ in the direction of anisotropy is 

{a^} = [c( a/? )]{e E ( a/? )} (4.11) 


where 
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iMoflx - iMafi EM E( a/?) 2 E(a/3) 2f E (a/3) E(a0) y T 
t e J“i e ll ,e 22 ’33 ’12 ’13 ’ 23 ' 


and 


[c( a/?) ] = 


,(a/3) 

'll 


c (°0) 

U 12 

Q( a P) 

^22 


c (°0) 

u 13 

Q (<*/?) 

u 23 

c( a ^) 

u 33 


0 

0 

0 


'44 


0 

0 

0 

0 


sym. 


,(a/?) 

'55 


0 

0 

0 

0 

0 

n( a P) 

U 66 


(4.12) 


Appropriate elastic constants are substituted in (4.12) to represent either fiber 
or matrix subcells. 

Combining (4.8), (4.9) and (4.10) the following relations for the average 
subcell stresses sj are obtained 

if 1 = cjf + c [f\i aP) + i aP) -^[f 

if = if *n + ifi a0) + ifi^-^Af 
if = if hi + ifi a0) + 4fM a/5) - ^ a fzf 


JaP) _ , A»ffh _ 2C ( <*%(<*« 

^12 “ U 44 ^7x7 J 2 ^44 12 


(4.13) 
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dvf 


S {(*P) _ Q {a0 ) r 3 ,{a0\ _ 2C {a0) L (a0) 

s 13 ” u 44 teT + H J 2U 44 l 13 


S (^) = C {°fi\A°0) + J< a 0) ]-2c( a ®lS a ® 

°23 fifi l™3 + j ^93 


"66 


66 23 


In (4.13) above, l| ^ represents the average total subcell plastic strains, and 
is the isotropic constituents’ elastic shear modulus. 

If classical rate— independent plasticity theory is selected to describe the 
nonlinear material behavior, then the plastic strain rates are given by 


i p M = A J><& 
ij a0 ij 


(4.14) 


which is the flow rule associated within the von Mises yield criterion, with = 

o- • - 8- -cr , , /3 representing the deviatoric stresses (£.. is the Kronecker delta) 
!J IJ IJ 

and A a p the flow rule function. The above leads to 



c(«0) 

b ij 


(4.15) 


where are the deviators of si^. Thus the plastic average strains li a ^ 
are determined by integrating the above flow rules. 

The conditions of traction continuity along the interfaces of the subcells 
in the first order theory result in 


g(l/2) _ g(2/?) 

5 2i “ 3 2i 
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(4.16) 


c( <*1) _ c( (&) 

5 3i _s 3i 

By imposing the continuity conditions previously stated, the variables 

(p( a ^ and can be eliminated. This leads to closed form solutions for the 

average subcell stress in (4.13), and the average stress in the cell o - in 

^ j ^ j 

(4.8). Without listing the details (see [64]), the results are given here: 
a ll = b ll £ ll + b 12 £ 22 + b 13 £ 33 ” H 11 

a 22 = b 22 £ ll + b 22 e 22 + b 23 £ 33 ~ H 22 

^33 = b 13 £ 1 1 + b 23 £ 22 + b 33 £ 33 “ H 33 

(4.17) 

a l2 = 2b 44 e 12 -H 12 
a 13 = 2b 55 e 13 -H 13 


a 23 " 2b 66 £ 23 H 23 


where the b- j are the effective elastic constants of the composite. The constants 
b- • and H- • are given in functional form as 
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b ij " f * C i j^’ v ofi h a’ V 

H ij - ^Sf 1 '• v v v (418) 


The resultant constitutive relations in (4.17) involve an overall elastic 
modulus matrix [B] consisting of 9 independent elastic constants b^, b 12 . b^g, 
b 22 , b 2 2 , bgg, b^, bg 5 , bgg (thus orthotropic). For square fibers and equal 
spacing the number of independent elastic constants reduce to six since b^ 2 = 
b^, b 22 = b^g, and b 44 = bgg. If a transversely isotropic material 
representation is desired, in which only five independent effective elastic 
constants are needed, then a transversely isotropic averaging technique 
described in [64] can be implemented. The overall constitutive relationship is 
summarized by 


{a} = [B]({?} - {i P }) 


(4.19a) 


where 

{i P } = [B]~ *{H} (4.19b) 

Equations (4.19) contain the average ply properties Bjj and the effects of 

_p 

applying plasticity at the subcell level, which are included in the term {e }. 
This is the full development presented by Aboudi. For the present work, only 
the transversely isotropic version of the B-j constants are of interest, as the 
plasticity is not implemented into the micromechanics theory. It is anticipated 
that this will be done in the future, and a systematic method for this 
implementation was developed by Arenburg [63]. In this study a macroscopic 
orthotropic plasticity development is included. 
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5. ORTHOTROPIC PLASTICITY FORMULATION 
5.1 Introduction 

In order to develop the capability of modeling orthotropic elastoplastic 
behavior from a macroscopic viewpoint, a yield function introduced by Sun et 
al. [74—77], which is quadratic in the stresses and employs the associative flow 
rule and isotropic hardening, was adopted. The plane stress radial return 
algorithm by Simo et al. [89,90] was modified to include this new yield function. 
This algorithm was chosen due to its accuracy, improved global convergence 
rates, and compatibility with the probabilistic finite element routines. In this 
chapter the basic governing plasticity equations will be stated, the radial return 
algorithm for their solution at each gauss point will be developed, and the 
solution of the finite element equilibrium equations along with the layerwise 
integration procedures will be discussed. 


5.2 Governing Equations 

According to Sun [76], the orthotropic yield function is given by 
f = ^ (a 11 a^ 1 + a 22 a 22 + a 33 a 33 + 2a 12 c7 ll cr 22 + 2a 13 a ll a 33 + 2a 23 a 22 a 33 


+ 2a 44 a 23 + 2a 55 <r 13 + 2a 66 a 12^ " Y 


(5.1) 


where Y is a state variable and a ■■ are the stresses in principal material 

I 

coordinates. The coefficients a^ are constants, which are determined from 
experimental data and control the amount of anisotropy in the plasticity. This 
yield criterion does not include the assumption of incompressibility of plastic 
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strains or that hydrostatic stresses result in no yielding or plastic deformation. 

The function also reduces to the von Mises criterion or the Hill yield function 

for orthotropic materials [67,68] with appropriate selection of the a- . values. 

^ J 

For example, for the von Mises function the a^ values are 



a 44 - a 55 “ a 66 “ 1 


The associative flow rule is employed, which allows the incremental 
plastic strains to be stated as 


de \i~ dl wr. ( 5 - 2 ) 

J iJ 

where d7 is an incremental plastic load parameter, and index notation has been 
used. It is the practice in plasticity to define the effective stress in such a way 
that it is related to the yield function and in a manner that it reduces to the 
stress in a uniaxial tension test. With this in mind, and following the 
development of Sun [76], the effective stress is defined as 

a = fit (5.3) 

The effective plastic strain must be related to the effective stress in such a 
way that plastic work considerations are consistent. Thus the quantity ade^ 
must represent the increment of specific plastic work of deformation 
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dW P = 


V'U- 


Equating the two quantities of plastic work, we have 


ade P = 


ij 



(5.4) 


Substitution of (5.1) into (5.2) and the result into (5.4) gives 


V'ij = 2 “ 7 

so that 

d jP = 2fd2 2 
a 


(5.5) 


It is also noted that Y becomes 


v 2-2 

Y = r 


(5.6) 


The shell element used here is degenerated from 3-D elasticity using the 
kinematics of the first order shear deformation theory. In other words, 
transverse normals remain straight and inextensible, i.e., = 0. 

Consequently, < 7^3 does not enter the strain energy of the shell. 

Introducing vector notation, the stress and strain tensors become 


{ a } - {cr l;L ^22 G \i ^13 ^ 23 ^ 


and 
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{ f } - { £ n e 22 2e 12 2e 13 2e 23^ ’ ~ { e ?l £ 22 2c ?9. 2e ?3 2e 9^} 


11 22 * c 12 13 23 J 


The components of the back stress q- • (included to model kinematic hardening), 

J 

and the relative stress tj^ = cr- • - q^, are expressed in the vector form as 


T T 

{q} = {Qh ^22 q 12 q 13 q 23^ ’ ^ = Wl ^22 ^12 ^13 ^23^ 


Thus the governing elastoplastic equations in the nonincrement al (time domain) 
vector form can be expressed as 


(5.7) 


M = { ( e j + < £ P) 

M = [C] {< e } 

or 

{fP} = 7 (associative flow rule) 

{,} = §h'{ £ P} 

f = 5 (a ll^ll + a 22^22 + 2a 12 T? ll 7? 22 + 2a 44^23 + 2a 55 , 13 + 2a 66 , 12^ = 

. 2 . _ 
a = 7 <7 

a = v/3T 


/ 

where 7 is the time derivative of the plastic load parameter, H is the kinematic 
hardening modulus, and the symbol a denotes the equivalent plastic strain. 
The matrix [C] is the elastic constitutive matrix, adjusted for the constraint 

a 33 = 0l 

When the flow rule in Equation (5.7)g is evaluated, the following result 
is obtained 
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0 p } = 7[P]{»} 


(5-8) 


where the matrix [P] is given by 

0 0 ' 

0 0 

0 0 

000 2a 55 0 

0 0 0 0 2a 44 

♦ 

The effective stress a can be expressed as 

» = J \ M T [P]« 

Using the matrix [P], the governing equations (5.7) can be recast as 


[P] = 


a ll a 12 0 

a 12 a 22 0 
n n 9a 


{<} = { e e ) + {(P) 

{a} = [C] { £ e } 

{iP} = 7[P]{»} 

{q} = 7jH" [P]{,} (5.9) 

P - 2 ( t ?} T [?]{*?} “j Y 2 (a) < 0 

i = 7 [j {»} T [P]{»}] 1/2 

where the parameter Y represents the hardening law in terms of the equivalent 
plastic strain a. 
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Loading and unloading conditions are stated in simple Kuhn— Tucker 
form [89] by requiring that 


For an elastic process, f < 0 and 7=0. For a plastic process f = 0 and 7 > 0. 
These two conditions are generally valid, whether in the loading or unloading 
state. 

5.3 Incremental Formulation 

Employing a backward Euler difference scheme to integrate equations 
(5.9) over time {t n ,t n+1 }, and letting 7 n+ ^ = 7 n+ ^At (the plastic load 

parameter) and f = V {v} [PK 7 ?}. the strain at t n+1 can be written in terms 
of the strain at t and the gradient of the incremental displacements, 


where "V" stands for the differential operator used in the definition of strains. 
A trial stress state is assumed by freezing plastic flow so that the entire step is 
purely elastic. The trial stress then becomes 


f < 0 , 7> 0 , 7f = 0 


(5.10) 




(5.11) 
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( c n+l^ - ^ e n^ + 1 'n+ll P ^’ , n+l^ 

{W = K> + W5 H 'l p ]<W 


(5.12) 


“n+l = °n + ' / 3 VhL+1 
Restating equation (5.12)^ as 




(5.13) 


substituting {Ae P +1 } = 7 n+1 [P]{>i n+1 ) and {i? n+1 } = {% +1 ) - <<ln++ 
and rearranging the terms, the following relations are obtained: 


where 


< W = i + h \ " 

1 + 3 n 'n+l 


[allCrM-lJli" 1 } 


PM 


icr 1 + 


r n+l 

1 + J H 7, 


[P) 


-1 


n+l 


(5.14) 


In order to perform the incremental updates required in equations (5.12) 
the plastic load parameter (Lagrange multiplier) 7 must be determined. It is 
found by enforcing the consistency condition at time t n+1> i.e., 
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(5.15) 


f2 ( 7 n+l) 5 \ ? n+l “ ^ Y(a n + Vn^n+l^ ~ 0 


The actual hardening functions used in the program are those recommended by 
Simo and Hughes [89] 

Y(a) = 0Ea 4- <7y 4- (K^ - ayMl — exp(— Aa)] (5.16) 

and by Sun [76], 


Y(a) = H 


a + 


A-. 1/A 


H 


(5.17) 


where H is the hardening modulus, ay the uniaxial yield stress, and A are 
other input parameters, and 

H(a) = (1 -/?)Ha 

Here /? denotes the fraction of kinematic and isotropic hardening desired, i.e. 
0 = 1 denotes purely isotropic hardening, and 0=0 denotes purely kinematic 
hardening. Equation (5.15) is solved at each gauss point in the structure for 7 
by a local Newton iteration procedure, as it is generally a nonlinear scalar 
equation. 

The global equilibrium is obtained by using Newton— Raphson iteration. 
This requires that the tangent moduli be known in the form 



Simo and Hughes developed tangent moduli, which are consistent with the 
integration algorithm previously discussed. For finite values of load step size, 
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they showed that the consistent elastoplastic tangent moduli preserved the 
quadratic rate of asymptotic convergence that is characteristic of Newton’s 
method. Use of the continuum tangent moduli derived independent of the 
algorithm loses this convergence rate. Differentiating the following incremental 
equations 

K+l> = [ C l({W-< ( »+l}> 

“n+1 = “n + ^1 Vl-l f n+l (5 19) 

{q n+1 } = {q n } + i'n + jH'[ p n,„ + i} 
and substituting equation ( 5 . 12)2 and = { cr n _j_ ^ } — we obtain 

< d VH> - ICHK+l) - d Vn™W - 7 n+1 [P]({^ n+1 > - {dq n+1 })] 

d “n+1 = /I < f n + l d Vn + WW ( 5 ' 2 °) 

0 / 

L H 

^ dq n+l^ = TZ ' ^ d7 n+l^n+l^ + 7 n+l^ da n+l^ 

1 + 3 M 7 n+l 

By regrouping equations (5.19)^ and ( 5 . 19)2 it can be shown that 
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< d<r n+l> = PI 


' , , d \+l 

( df n+l> _ ^TT 


1 + T H X 


[ p l<W 


(5.21) 


n+1 


{d, "« } = 7TTi^ 


< d< 'n+l>"l H ^n+l^n+l* 


5 “ 'n+1 


Differentiation of the consistency condition (5.15) and use of the definition of 
f ^ results in the expression 



7 n+l^n+l} l P K d Vl-l 



d7 n+l 


= 0 


(5.22) 


Substituting equations (5.20)^ and (5.21) into equation (5.22) and solving it for 
d 7 n _|_i we obtain 


¥W T l p il ; H d W 

( 1+ ^+iH”n + ii 1 '[ p n s n p Hi„ + i} 


(5.23) 


where 


*1 - 1 + 5 H 7 n+l 


°2 ~ 1 “3 Y 7 n+l 
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_ 2 « [Y'^ + H'g 2 ] 

n+1 * T 2 n+1 {W T [ P ]l 2 )[ P KVM> 

Substituting (5.23) into (5.20^ and using the definition of [S], the consistent 
elastoplastic tangent moduli can be expressed as 


5.4 Computer Implementation 

For the general case of orthotropic plasticity, the constitutive matrix [C] 
is given by 


'11 


D 


v \ 2 E 22 
D 


0 0 0 


U 12 E 22 E 22 
D D 


0 0 0 


[C] = 


G 12 0 


0 0 
0 0 


0 KG^2 0 

0 0 KG 23_| 


(5.25) 
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where D = 1 — ^ 2 U 21 anc ^ ^ ' s t ^ ie s ^ ear correction factor taken to be 5/6. If 
yielding occurs, it becomes necessary to solve the scalar consistency condition of 
equation (5.15) for the load parameter 7 n _^- A Newton iteration is used here 
due to the general nonlinearity of (5.15). Thus, if the variable F is used to 
represent the consistency condition, then the Newton iteration is performed by 
employing 


Vi 



(5.26) 


at each yielded gauss point. Using the chain rule of differentiation, the 
necessary derivatives of (5.15) can be obtained in a straightforward manner for 
the general case of orthotropic plasticity. 

For the special case of isotropic (von Mises) plasticity, Simo and Hughes 
[89] have developed a simplified form of the consistency condition that is a 
direct function of the load parameter 7 j.. This extra work results in an 

effective reduction of the number of computations required in the Newton 
iteration. Recall that for isotropic elasticity, the matrix [C] is modified by 
= E 2 2 = E, 1^2 = ^21 = u > anc ^ G>i2 = = G 2 3 = G, and the matrix [P] is 

modified due to the von Mises values for the coefficients a-j previously 
mentioned. Thus the isotropic matrices [C] and [P] can be simultaneously 
diagonalized using the following spectral decomposition 


[P] = [Q][A p ][Q] T and [CJ = [Q][A c ][Q] T 


(5.27) 
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where 


[ 1-1 0 0 0 ‘ 

, 110 0 0 

[Q] = — 0 0 ^ 0 0 

0 0 0 0 
0 0 0 0 y/7 



'1/3 0 0 0 O' 

0 10 0 0 
0 0 2 0 0 

0 0 0 2 0 

0 0 0 0 2 



E 

l—i/ 

0 

0 

0 

0 


0 0 0 0 

2G 0 0 0 

0 G 0 0 

0 0 GK 0 

0 0 0 GK 


Using the matrix [Q], the following transformation is introduced 


{£} = [Q] T to) = 


’ton + ' 

H 11 + 7 7 22 )/v2 
77 12 
77 13 
77 23 


The mapped trial state in terms of the elastic trial state defined in equations 
(5.12 ) 2 can also be transformed as 


60 


(5.28) 


r ,t r ial 


} = [Q] T {x 


t r ial 
n+1 


} 


Substituting relations (5.27) in (5.14) the diagonalized update equations can be 
expressed as 

^n+1^ = 




0 0 0 0 


^+2G 7 


0 0 


0 0 


0 0 


^+2GK 7 


0 0 0 


^+2GK"7 


f t r ial-i 

K+l } 


(5.29) 


Transforming the consistency equation into {£} variables, and utilizing equation 
(5.29), we have 
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*cw-» 


ft r ial . triak2 /trial triak2 , t r i a K2 

l^ll +^22 ) , 1^22 _7? 11 ) .^12 ) 


■° ^l + 3(l-i *)} 


+ 


( 0 , +2GiY 


+ 


(0!+2G 7 y 


ft i i a K2 / t r i a 1 \2 ■ 
, <*13 ) , ^2 3 ) 

(^+2GK 7 ) 2 (^+2GK 7 ) 2 


-^Y 2 (a 


n+1 


) = 0 


(5.30) 


This is the form of the consistency equation that is used in the scalar Newton 
iteration at each gauss point to solve for 7fl+1 for the case of isotropic von 
Mises plasticity. 

The radial return algorithm for the constrained = 0 state of stress 
can be summarized as follows: 

1. Update the total strain using (5.11) and compute the trial elastic 
stresses from (5.12)^ 2 - 

2. Compute f^J d from (5.15) or (5.30). Exit if f^j 31 < 0; 
otherwise solve consistency equation (5.15 or 5.30) for 7r+1 using 
Newton iteration. 

3. Compute the algorithmic moduli matrix [E] from equation 
(5.14) 2 . 

4. Update relative stress, back stress, actual stress, plastic strain, 

and equivalent strain from (5.14)p (5.12)^, {o^^} — + 

{q n+1 }, (5-12) 3 , and (5.12) 5 , respectively. 

5. Compute consistent elastoplastic tangent moduli from (5.24). 
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6. Update e,^ using, if desired, 


13 


23 


33 n+l ^n+l + ^22 ^n+l- 1 


P , P 

11 ,, + e 22 , - 
n+1 n+1 


5.5 Layerwise Integration Procedure 

Recall equations (3.19), (3.20) and (3.22): 


i k l i - . 

v [B L ] T [Cl[B L ]dV 

(3.19) 

t - 1 

II 

e 

v t B NL^ S ^ B NL^ dV 

(3.20) 

{F} = 

J v [B L ] T <S}dV 

(3.22) 


From equations (3.19) and (3.36) [Kj can be written explicitly in terms of the 
thickness coordinate ( as 

[KJ = [ ([SD1] + C[SD2] + C 2 [SD3]) T [C , ]([SD1] + C[SD2] + £ 2 [SD3])dV 

J y 

(5,31) 

Gauss quadrature is used to perform the numerical integration. In Gauss 
quadrature an integral over the interval [—1,1] is replaced by the weighted sum 
of function values at selected points and weights: 
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1 


-1 


N 

f(o d £= s w i f (y 

i=l 1 1 


(5.32) 


where Wj are weights and (■ are the base points. Equation (5.31) is evaluated 
using gauss quadrature: 


P 3 3 

[K T ]= £ £ £ 

k=l £=1 77=1 


2 

£ 

C=1 


1 b l 1 l c !(k)l B Ll w f w ^ w c°l J l 



(5.33) 


/ 

where [BjJ denotes the expression in the parenthesis of equation (5.31), and 
P = total number of layers in the laminate 
h^ = the thickness of the kth layer 
h = total thickness of the shell 
0 1 J | = determinant of the Jacobian matrix 
For elastic structures, the problem is reduced to a 2— D one by performing 
explicit integration through the thickness. The Jacobian matrix, in general, is a 
function of £, 77 , and the ( terms in the Jacobian matrix may be neglected if 
the thickness to radius of the shell ratio is small. If the Jacobian °|J| is 
independent of (, explicit integration can be used. In addition, we neglect the 
( terms in [BjJ. With these assumptions, [KjJ in equation (5.31) can be 
written as 
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[K l ] = E E [[SD1] T [A1][SD1] + [SD2] T [A2][SD1] 


+ [SD1] T [A2][SD2] + [SD2] T [A3][SD2]j 0 1 J | 


w> w 


£ W ^ 


(5.34) 


where 


[Al] = 


P 

E 

k=l 


A+i 


[C']dC 


P f C k+l , 
[A2] = E C[C }d( 

k =i J c k 


[A3] = 


P A+l 

E 

k = iJ c k 


C 2 [c']dC 


(5.35) 


It is important to note that the upper limit on the number of quadrature points 

is set to the value 3, which implies full integration. Reduced integration implies 

/ 

a value of 2. For elastic structures, [C ] is constant within each layer and thus 

(5.35) can be integrated explicitly. For elastic-plastic structures, some 

/ 

locations through the thickness become plastic, so that [C ] must be replaced by 
[C e P], the elastoplastic tangent moduli. Thus, when plasticity is included, 
numerical integration must be used through the thickness as well. Gauss 
quadrature will also be used in the thickness direction with the thickness 


65 


coordinate given as 


C - C k + ir + z i^ 


(5.36) 


Here is the base point, and — 1 < z < 1 for each layer. Substituting equation 
(5.36) into (5.35), then [Al], [A2], and [A3] take the form 


p N h, 

[Al] = £ E [C] Yi w i 
k=l i=l 1 


p N h, _ h, 

t A 2 l= s . . s . ^k + f~ ^ + z iW c l r w i 

k=l 1=1 


(5.37) 


p N h, 2 _ h, 

[A3]= E E (C k + ir( 1 + z i^ ^H“ w i 
k=l 1=1 


where 

[C] = [C ] for an elastic gauss point, or [C ] for a plastic gauss point 
N = number of gauss points per layer 

Next we consider the matrix [K^jJ from equations (3.20) and (3.36), 
which can be rewritten as 


[DHl] 


p r 
s 

C=1 J 


>k+l 


[S ]dC 


[DHl] 
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+ [DHl] 


S 

k=l 


A+l 


C[S ]dC 


[DH2] 


+ [DH2] 


P 

S 

k=l 


A+l 


CIS ]d( 


[DHl] 


+ [DH2] 


' P 
S 

k=l 


A+l 

^k 


CIS ]dC 


[DH2] 


°| J | W w 

s r 's 


(5.38) 


The stresses [S ] are calculated from either the elastic constitutive relations or 
the radial return algorithm when the gauss point yields. The integrals 


P 

[N] = E 
k=l 


A+l , 


Cl 


P A+l , 

[S ]d( , [M] = E C[S ]dC, and 

k=l J Ci 


[M2] = S 
k=l 


5 k+l 


CIS ]dc 


(5.39) 


must be evaluated using gauss quadrature as in (5.37) when plasticity occurs. 
The load vector {F} from equations (3.22) and (3.36) can be written as 


{F} = J j[SDl] T + C[SD2] T + ( 2 [SD3] T J{S }dV (5.40) 
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o 

Regrouping the terms and neglecting ( , we can write 




(5.41) 


The thickness direction integrals are evaluated numerically when plasticity 
occurs: 


p A+i , 

{N} = £ {S }d< 

k =1 Jc k 


P c^k+1 , 

{M} = £ <{S }d< (5.42) 

k= lJ C k 


This numerical integration through the thickness for plasticity greatly increases 
the computational expense. 


5.6 Solution of Equilibrium Equations with Plasticity 

p 

At each load step, the plastic strain {e n } and the back stress {q n } are 
"frozen" for the next (n+1) load step and are not updated until convergence is 
achieved. This is important so that inaccurate iterative plastic strain and back 
stress increments are not used. 
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The Newton Raphson method described for the geometric nonlinear case 
is also used here to solve for global equilibrium. When a gauss point yields, the 
elastoplastic tangent moduli are substituted. Two convergence criteria have 
been installed: 1) a displacement norm, and 2) a force residual norm. 
Convergence is achieved if the specified norm is within a preassigned tolerance 
e, usually selected to be .001 or less. The displacement norm is defined as 


'NEQ 

E 

J = 1 


(A 


0 )_ 


/j_i\ o 

A* 1 1 >YI E 
J j = l 


< A i°> 



e 


where NEQ is the total number of equations (or degrees of freedom) in the 
model. The force residual norm is given by 


'NEQ 

E 

.1 = 1 


rW _ FV 
J J 


(i) 


0 NEQ 9 

2 / £ (RO 2 

j = l J . 


1/2 


< e 


where R- is the jth component of the external load vector {R}, F. is the jth 
j j 

component of the internal force vector {F}, and the superscript (i) denotes the 
iteration number. 
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6. SECOND-MOMENT PROBABILISTIC FINITE ELEMENT METHOD 

6.1 Introduction 

In this chapter, the second— moment probabilistic theory for linear, 
geometrical nonlinear, and material nonlinear analysis are developed for time 
independent behavior. The selection of random variables are discussed, along 
with methods for computing the necessary derivatives. Computational savings 
techniques are presented, and various correlation assumptions discussed. The 
theoretical developments given are as applied to the degenerated 3— D shell 
element discussed in Chapter 3. 

6.2 Second— Moment Probabilistic Method 

In Chapter 3, the finite element incremental equations were cast as 

[K]{A} = {R} (6.1) 

which is applicable to both linear and nonlinear analysis. Following the 
development of Liu et al. [33 — 40], equation (6.1) can be rewritten as 

{F({A},{b})} = {R({b})} (6.2) 

where {F} is the internal force vector, {A} is the displacement vector, {R} is 
the external force, and {b} is a discretized vector of the random function b(x), 
where x is a spatial coordinate {x}. As in typical finite element analysis, the 
random function b(x) is expanded using shape functions ip-^x) 

b(x) = S V>.(x) b- (6.3) 

i=l 1 1 
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where bj are the nodal values of b(x). Generally the random quantity b can be 
a material property, geometric dimension, or a load. 

The second— moment probabilistic method can be mathematically viewed 
as a perturbation method, in which only up to second— order terms are retained. 
Following the typical perturbation procedure, the matrices in equation (6.2) are 
expanded in Taylor series about the mean value of the random quantity of 
interest b. After substituting the Taylor series into (6.2) and equating similar 
order terms, one obtains equations involving various order derivatives of the 
displacements with respect to b. These displacement derivatives can be used to 
determine the statistical mean and covariances of both the displacement and 
element stress. 

In order to describe the probabilistic distribution of the input random 
function b(x), two statistical moments are required as input. In addition, if 
spatial correlation is used, this functional dependence must be described as well. 
Using the notation in [33], the expectation b(x), denoted by E[b(x)], coefficient 
of variation a, and autocorrelation coefficient function A(b(xj), b(xj)) must be 
defined. These quantities are known inputs to the model. Typical 
representations for the autocorrelation coefficient function are of the form, 

A(b(x i ), b(xj)) = exp(— | - x • | / A) (6.4) 

where A is the known correlation length of the random field. The covariance 
cov(b-,bj) is given by 

cov(bj,bj) = [var(b(x i ))var(b(xj))] 1 / 2 A(b(x i ),b(Xj)) (6.5a) 


where 
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varCbCxj)) = a 2 E[b(xj)] 2 (6.5b) 

The definitions of autocorrelation and the actual functions used are discussed in 
more detail in the sequel. 

The vectors {A}, {R}, and {F} are expanded about the mean of the 
random function b using the Taylor series: 


n 


n 


{£} = {£}+ E {£}.db i + i E {A}. . dbidb- 

i=l D i 1 ^i,j=i D i D j 1 J 


{R} = {R}+ E {R}.db.+i E {Rk.db.dbj 
i = l °i 1 %j = l D i D j 1 J 


n 


{F} = {F} + S {F} b + [K T ]{A} b 
i=l ^ i i 


db: 


n 


+ E 

ij=l 


\ { p h).b. + 5 [K T ]{A} b b + l gT ]b.^ A ^b. 

1 J 1 J li 


dbjdbj 


( 6 . 6 ) 

(6.7) 


(6.8) 


In the above equations the following notation is used: 

b(x) = E[b(x)] — for mean value or expected value of b 

dbj = e(b- — b.) — first order variation of b- about bj 

§b = ^ — partial derivative of a function g with respect to bj and 

evaluated at 6- 
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Also [K T ] 


is defined to be the tangent stiffness matrix, 



Use of the above approximations in equation (6.2) results 
following perturbation equations: 

Zeroth— order equation 

{F} = {R} 

First-order equations 

[R T l{A} b ={fi} b -{F} b , i = l n 

i i i 


Second— order equations 

[k t ]{a 2 } = {jy 

where 

i A 2> = 2. E ^ A >b.b cov ( b i> b j) 

i,J=l i J J 

and 


{R 2 >= e 

l,j=l 



- [K T ] b {A} b Icovtbpb.) 
1 j J 


( 6 . 9 ) 
in the 


( 6 . 10 ) 


( 6 . 11 ) 


(6.12a) 


(6.12b) 


(6.12c) 
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Once {A}, {A}^ , and {A2} are obtained by solving equations 

(6.10)— (6.12), the mean and autocovariance matrices for the nodal 
displacements can be determined. These are formally defined as 


E[{4}] = 


{A({b})}f({b})d{b} 

J 0D 


(6.13) 


and 


cov(A r ,A s )=f° (A r - A r )(A s - A s )f({b})d{b} (6.14) 

** — 0D 

respectively. Here f denotes the joint probability density function, A 1 is the ith 
degree of freedom of {A}, and {b} is the random variable vector. By 
substituting the Taylor series expansion of {A} from equation (6.6) into (6.13), 
the second— order estimate of the mean value of {A} is obtained (see [23]), 

E[{A}] = {A}+1{ E (A). . covfbj.b,)) (6.15) 

i,j = l i j J 

Similarly, the first— order cov(A 1 ,A''), which is consistent with a second-order 
analysis, is given by 


cov(A r ,A s ) = S A' A? cov(b.,b.) 

i,j=l i j J 


(6.16) 
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Recall that the stress vector in an element is given by 

M = [C]{ f } = [C] [B]{ A} (6.17) 

The functional dependence of {a} on {b} is through the constitutive matrix, 

[C] = [C(x,b)] 

and [B] is the linear or nonlinear strain— displacement matrix. The Taylor series 
expansion of [C] is 


[C] = [C] + s 
1=1 


n 


{C} b dbj + £ 


i,j=l 


{C} b .b db i db j 

i j J 


(6.18) 


Substituting (6.18) and (6.6) into the definition of E[a] (similar to equation 
(6.13)), the second— order mean stress can be written as 


E[{ cr}] = [C]E[{e}] + 


n 

S 

i.j = l 


[C] b .([B]{A}) b . 




cov(b i ,bj) 


(6.19a) 


where 


E[{«}J = IB]{A} + \ 


n 

2 

.iJ =1 


([B]{A}) b b cov(b i ,b.) 

i j J . 


The autocovariance of stress can be expressed as 


(6.19b) 
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( 6 . 20 ) 


n 

cov({a e },{<T f }) = _ V { <7 e }b.{ £T f}bj COV ( b i> b j) 

where {a g } represents the stress vector for element e. The vector {a is 

u i 

obtained from (6.17) as 

{ ff e>bj = I C e) bi l B eH A e) + t«y([B e ]{A e » b . < 6 - 21 > 

Thus (6.20) can be evaluated using (6.21). 

6.3 Composite Random Variables 

Sources of randomness can be material properties, geometric dimensions, 
or loads. For the present study, the only geometric dimensions selected as 
random variables are the ply thickness and ply angle. The loading is considered 
to be deterministic throughout this study. All material properties are treated 
as random variables, either from the point of view of the ply level or the micro 
level (when a micromechanics constitutive theory is incorporated). At the ply 
level, material variables could be any of the engineering material properties 

Eff, E 22 , ^ 2 ’ ^12’ ^13’ ^23‘ ^ eca ^ ^ at l ^ ese properties along with the ply 

/ 

angle 6 define the coefficients of the constitutive matrix [C ] in equations (3.34) 
and (3.35). At the micromechanics level, the material variables could be the 
fiber and matrix properties: E^, E^i ^fi2> V {12* u f23’ ^m’ ^m’ w ^ere 

"f" denotes fiber property, "m" denotes matrix property, and FVR is the fiber 
volume ratio. These micro— variables are used in the Aboudi micromechanics 
equations to determine the engineering (ply level) material properties. In this 
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section, detailed discussions of the technical assumptions and subtleties of 
incorporating these random variables are presented. 

6.3.1 Ply Thickness 

The total thickness of a laminate is determined by summing the 
individual ply thicknesses, which for a composite made of one material type, are 
assumed to be a constant for every ply. Here it is assumed that the thickness of 
all layers fluctuate, but that the total thickness of the laminate (shell) remains 
unchanged. This is the assumption made by Nakagiri et al. [26] for eigenvalue 
analysis of composite plates. In [26] only ply thickness and orientation angle 
were selected as random variables. Since the total thickness is assumed to be 
constant, the ( coordinates of the upper and lower surfaces of the shell are 
deterministic. The through— the— thickness integration, as described in 
equations (5.35), involves the ply thickness coordinates as follows: 


p A+i , 


[Al] = S 


[C ]dC 




<(c']d< 


[A2] = S 




<V]dC 


[A3] = E 



For elastic structures, these integrations are carried out explicitly, 
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( A1 ] = k M c '] k «k + i-Ck> 

l A2| = j jM c \ ( <k+i-<k) («- 22 ) 

l«l-j k Mc'] k (<k + i-<k) 


where is the lower coordinate of the kth layer. For elastic— plastic structures, 
the integrals must be numerically evaluated. For a two— layer composite, for 
example, the previous assumption results in only the upper coordinate ^ 
layer 1 being random. In general, Cfc+i is the random variable for the kth 
layer, unless k is equal to P, the last layer. 

The random variable is also interpolated using the finite element 
interpolation functions V'j(x): 


n 


^k+1 “ ^¥ x )^k+l)i 


(6.23) 


Substituting (6.23) into (6.22), the perturbation derivatives {F}^ required in 

equations (6.11) can be evaluated for the thickness random variable. For the 
kth layer, assuming again an elastic constitutive law, the results are 
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(6 24) 


<S=([ C V[C'WlM 

4^ ) := ( [c' 1 k -[c'] k+I ) c k+1 ^ i 

4^: = ([c'] k -[c'] k+1 )O i 


Since the matrices [Al], [A2], and [A3] contain the only dependence on the 
random variable ^ equations (6.24) are used directly to assemble the 
resulting {F}^ . 


6.3.2 Ply Angle 

The only dependence on the ply angle 6 exists in the constitutive matrix 

/ 

[C ] (see equation (3.54)). Therefore, in order to evaluate {F}^ , where b is the 

ply angle, it is necessary only to determine the derivative matrix 1 for each 

layer. This differentiation is done explicitly and is straight-forward, so the 
details are not presented here. 


6.3.3 Ply— Level Material Properties 

The only dependence on the ply— level engineering constants E^, E 22 , 

/ 

^12’ ^12’ ^13’ ^23 * s * n e( l uat * ons (3.34) and (3.35). The derivatives J , 

u i 
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where b is any of these properties for a particular layer, are explicitly 
differentiated here as well. 


6.3.4 Micro-Level Material Properties 

The micro— level random variables E^, E^) ^fi2> v {\2' ^123 > 
and FVR are specified inputs to the Aboudi micromechanics equations with the 
output being the ply— level engineering constants. These engineering constants 
are used to calculate the [Al], [A2], and [A3] matrix stiffnesses in equation 
(5.35). It should be noted that when the micromechanics model is used, only 
linear or geometrical nonlinear behavior is allowed. Thus, in order to evaluate 
the {F}^ matrices needed in equations (6.11), equations (6.22) must be 

differentiated in terms of the micro— level random variables. Referring to 
equation (6.22)p the chain rule of differentiation is performed 


d[ All _ d[ All dE l 1 , df All dE 2 2 , d[Al] dl/ 12 

3fVR ” dE^dFVR + dE^ 3FVE + 


, d [ All dG 1 2 , d f All dG 1 3 , d[All dG 23 

+ dFVK + dFVK + dTVE 


(6.25) 


where the selected micro random variable is FVR. Note that the terms d j^^ , 

d j^^ etc. are already known explicitly from the discussion in section 6.3.3, so 
dE. . dE 2 2 

that only the terms dFVR’ be computed. 
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The relationship between and FVR is given in terms of the 

micro mechanics equations. Hence, these equations must be differentiated in 
dE 

order to determine ^ pyp . Since the Aboudi model is rather complex, the finite 
difference technique is used. The derivatives are calculated as follows: 


Delta is a small number, usually set to .05 or less. In this way the perturbation 
derivatives {F}^ for micro random variables are successfully calculated for 

linear elastic or geometrical nonlinear elastic structures. 

6.4 Linear and Geometric Nonlinear Elastic Problems 

In this problem group, all of the necessary perturbation derivatives at 
the ply level are determined by exact differentiation of the stiffness matrices in 
equations (6.22), due to the elastic behavior. Finite difference derivatives in 
(6.26) are still used for micro— level random variables. For linear elastic 

problems, the vector {F} in equation (6.10) is represented by [K]{A} (as in 
(3.23)) and the derivative matrix {F}^ by [K]^ {A}. For geometric nonlinear 

elastic problems, {F} is given by: 


FVR'*’ = (1 + Delta)FVR 
FVR - = (1 - Delta) FVR 


( 6 . 26 ) 



where E^ 


is evaluated using FVR’*’ and E^ is evaluated using FVR , and 


81 


{F} = [ [B L ] T {S}dV 
j y 


Therefore {F}^ becomes 


<F} b . = f [B L ] T {S} b dV (6.27) 

i J V i 

The solution procedure involves a consecutive solution of equations 
(6.10) through (6.12). After the deterministic zeroth— order equation (6.10) is 
solved, either once for the linear case or iteratively at each load step in the 
nonlinear case, the generalized displacement vector {A} is used to perform the 
perturbation solutions in equations (6.11)— (6.12). In (6.11), there are as many 
solutions for {A}^ as there are number of nodes in the model, and in (6.12) one 

solution for ^2)- In addition, the computations in equations (6.11) and (6.12) 
must be performed for each layer in the model, as a particular random function 
is assumed to be independent from layer to layer. If there are n nodes and P 
layers in the model, (n + 1)P more matrix solutions are required at each load 
step for each random variable. This is not as expensive as it seems, because the 

_ rr 

stiffness matrix [K ] is inverted once and used in equations (6.11) — (6.12). 

The next step is to determine the mean and variance of the response. 
This is done using equations (6.15)— (6.16) for displacement and (6.19)— (6.20) 
for stress. Note that the derivatives of the term [B]{A} are required in 
equation (6.19) and for the stress derivatives in equation (6.21). For linear 


82 


structures this becomes 


([B L ]{A}) bi = [B L ] b; {A> + [B L l{A} bi (6.28) 

and for geometric nonlinear elastic problems it is 

«B NL ]{£}) bi = [B NL ] b .{A} + [B NL ]{A} b . (6.29) 

which implies that the left hand side of the equations are functions of both {A} 
and {A} b . 


6.5 Material Nonlinear Problems 

As in the geometric nonlinear problem, the residual vector is given by 
{F}. Therefore, the derivatives are the same as in equation (6.27). Recalling 
equation (6.2), {F} is functionally represented in the form 

{F} = <F({A},{b})} (6.30) 

Temporarily dropping the vector braces, and differentiating F with respect to b 


dF _ 5F , <9F 5A 

cJF 15E cHT 


(6.31) 


The first term on the right hand side involves the explicit derivative of F with 
respect to b. This is the result required in equation (6.27), in which the 
derivative of F is expressed in terms of explicit derivatives of the stress vector 
S. As stated in the last section, all these derivatives can be expressed exactly 
when elastic behavior exists. However, when material nonlinear behavior is 
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present, finite difference derivatives are used once again due to the complex 
nonlinear relationship between stress and strain. Since orthotropic plasticity is 
present in certain layers, the radial return algorithm discussed previously is 
considered to be very good for this purpose [38,40]. This is true since at a 
particular load step, the plastic strain and effective plastic strain are "frozen" 
and only updated after equilibrium conversion is achieved. This update can be 
done after the finite difference calculations are made, so that the true effect of 
perturbing the random variable about its mean is measured. Other advantages 
of the radial return method are the increased accuracy involved in the stress 
recovery routine and the algorithmic compatible tangent moduli which results 
in a more accurate tangent stiffness matrix. This accuracy is important in 
computing the perturbation derivatives (e.g., {A}^ in equations (6.9) — (6. 12)) . 

The evaluation of {F}^ involves the explicit differentiation of {S}, 

which is achieved using the finite difference formulas in (6.26), as recommended 
in [38,40]. Define 

bt = (1 + Delta)b 
bj = (1 — Delta)b 

S + = S(bt) and S _ = S(b[) (6.32) 

so that 

dS M S+-S" • _ , „ 

W = SlUeTtiJbj 1 “ 1 ”’ n 

It should be noticed in equations (6.32) that the perturbations on b are made at 
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a particular nodal location i to evaluate while for all other locations the 

derivative is zero. In this manner the applicable random variable derivatives 
can be evaluated. These include the ply level stiffnesses (E^, E 22 , G^, 

G 23 , G 12 ), ply thickness, and any plastic parameters desired. The plastic 
random variables selected are the uniaxial yield stress (Jy an d the hardening 
modulus H. 

In order to evaluate the autocovariance of stress, equation (6.20) is used 
as usual in conjunction with the total derivatives of S: 


or 


dS _ dS , dS bA 
db bb bC bb7 


dS _ dS ,dSde 
db7 bb be bb 


(6.33) 


The first term on the right hand side of ( 6 . 33)2 is the explicit derivative already 

or 

evaluated using finite differences. The second term is simply the tangent 
moduli matrix given by equation (5.18), and (bringing back the vector braces 
notation) {e} is represented from (6.17) as 

M = [B]{ A) (6.34) 

These derivatives are the same as expressed in equations (6.28) and (6.29) for 
geometrically linear and nonlinear structures, respectively. Substituting the 
results into (6.33), we obtain 
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( 6 . 35 ) 


{S} b . = if + [C eP ]([B]{A}) bi 

which is used to evaluate the autocovariance in equation (6.20). 

When combined geometric and material nonlinear behavior exists, the 
previously mentioned methods of differentiation are used for the appropriate 
plies. That is, exact differentiation is used for general elastic layers, and finite 
difference for elastic— plastic layers. The finite difference method is also used for 
the micro— level random variables, which are not applicable to the 
elastic— plastic layers. 

Up to this point the second— order equations have been developed. It 
turns out that for layered composite materials, independent random fields have 
been modeled for each layer. This greatly increases the cost of both assembly of 
the perturbation derivatives and of the solution of the first and second— order 
equations. Due to the immense assembly and storage requirements as well as 
computational effort required for the second-order equations when multiple 
layers exist, the analysis has been limited to the first-order equations. It can 
be observed that the second— order term results only in a small correction to the 
mean displacement and stress response, and it is assumed negligible when 
weighed against the cost of obtaining it. 
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6.6 Spatial Correlation 

The probabilistic finite element procedure developed herein has the 
ability to model the correlation involved in spatial fields. Consider the random 
variable vector {b(x)}, which can be a general random field, e.g., a family of 
related random variables such as Young’s modulus which can vary in the spatial 
coordinate x. The autocorrelation coefficient function of the discrete random 
field {b(x)} is thus defined similar to equation (6.14) as 


A(b(x i ),b(Xj)) = 


^ (b. -b i )(b j -6 j )f({b})d{b} 
[var(b(x i )) var^x^))] 1 / 2 


(6.36) 

The term "auto" here refers to the fact that we are dealing with the same 
random field {b}, and are concerned with its correlation in space. In order to 
determine the autocovariance of the response such as displacement and stress, 
A(bj,bj) must be a known input characteristic of the input random field b. A 
typical representation for the autocorrelation was given in equation (6.3), but 
this case was for isotropic spatial correlation. Extending this to an orthotropic 
ply and referring to Figures 6.1 and 6.2, we assume an orthotropic 
autocorrelation coefficient function in the following form 


A(b i ,b.) = exp 


/ / / / 



r 

1 

V 

2 

exp 


[V’jl 

2 


A,' 



A ' 



z J 

. 



V J 

_ 


ij = l-,n 


(6.37) 

In Figure 6.1, £ and 77 are the curvilinear coordinates of the shell which are 

/ / 

aligned with the body coordinates, £ and 77 are the curvilinear coordinates in 
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Figure 6.1 Local coordinate systems for body and material coordinates. 
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Figure 6.2 Body and material coordinates in the shell curvilinear plane. 
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the principal material coordinate system, aligned with the fibers, and A^,, A^, 

are input correlation lengths. Transforming £ and 77 to the principal material 

/ / 

coordinates £ and 77 , we obtain 




cos 6 

sin $- 

;«i-y 

l(*?i -*?jX 


—sin 6 

cos 6 

- *?j). 


(6.38) 


The formulation for the autocorrelation allows ply random variables to be 
correlated in the plane of the shell aligned with the material coordinates. For 
the case of material properties, we assume that they will have correlation trends 
which are parallel and perpendicular to the fiber directions. 

An underlying assumption in equation (6.37) is that A(bj,bj) is not a 
function of the shell thickness coordinate (. This leads to independent random 
fields from ply to ply (uncorrelated). This assumption is not a requirement, 
and correlation could be assumed in the ( direction, but for all of the examples 
given in this work all layers have been chosen to be independent. 

If the form of equation (6.37) is studied in relation to the correlation 
lengths A^' and A one discovers that for small A^' and A^' (say 1% of the 
length of the shell) that very little correlation A(b-,bj) exists from node i to j 
and corresponds to almost independent random variables. On the other hand, 
for large A^' and A^' (say greater than four times the length of the shell) the 
correlation becomes almost constant from node to node. As will be shown in 
the next chapter, this is equivalent to assuming a "uniform variance" field in 
which the fluctuation due to each random field is uniform over the entire shell 
layer. Basically this random field is fully correlated, and is equivalent to 
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assuming a single random variable for that ply. This is the same assumption 
made by Nakagiri [26], and it considerably simplifies the computations. 
Whereas previously the first-order equations (6.11) must be solved for all n 
nodes and for each of the P independent layers, the uniform variance 
assumption requires (6.11) to only be solved for each layer. This greatly 
reduces the number of computations required in (6.11), (6.16) and (6.20) as the 
loops extend only to P layers as opposed to P layers times n nodes. Thus if the 
random field is highly correlated so that the variance is essentially the same 
from node to node, then the random field can be dropped and a uniform 
variance field assumed to reduce the cost of the analysis. 

6.7 Computational Saving Techniques 

Three methods are utilized in this work to reduce the computational 
effort involved when a nonuniform correlated field is assumed. Liu et al. 
[33-40] discussed a method for diagonalizing the covariance matrix, and a 
technique for computing only the kth component of the displacement 
derivatives, which they called the adjoint method. Both of these are 
incorporated here, as well as a way to reduce the assembly time for stiffness and 
residual force derivatives by using the chain rule of differentiation. All three of 
these techniques will be discussed in this section. 

6.7.1 Diagonalization of Covariance Matrix 

In order to compute the covariance of the displacement and stress 
response, the multiplications involved in the double summations on i and j in 
equations (6.16) and (6.20) must be performed. As mentioned in [33], the 
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number of multiplications is proportional to n(n + l)/2, and is expensive. 
However, if the double loops on i and j could be reduced to a single loop, then 
substantial savings would result (the number of multiplications is reduced to n). 
This is accomplished by transforming the covariance matrix cov(bpbj) to a 
diagonal variance matrix var(cj,Cj) so that the following conditions hold 
var(CpCj) = 0 for i * j 

and 


var(Cj,Cj) = var(Cj) for i = j (6.39) 

The transformation is performed by solving the following eigenproblem 

[G]M = MIA] (6.40) 

where [G] and [A] represent cov(bj,b-) and var(cj,c-), respectively, [<p] is the 

J J 

matrix of eigenvectors normalized according to 

= ra 

[A] = W1 T [GM (6.41) 

and [I] is the identity matrix. The random variable vector {b} is transformed 
to the new random variable vector {c} according to 

{C} = M T {b} 
or 

{b} = M{c} (6.42) 

The perturbation equations (6.10) and (6.11) to be solved take the form (only 
dealing with first-order) 

Zeroth-order 

{F} = {&} (6.43) 
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First-order 


[K ]{A} C ={R> C -{F} c , i = l,..n 

i i i 

and the covariance solutions in equations (6.16) and (6.20) become 


(6.44) 


cov(A r ,A s ) = E A* A® var(c-) 

j=l j j J 

(6.45) 

n 

cov({a r },{<7 s }) = E {a r } c {a g } c var(cj) 
J ^ J J 

(6.46) 


Here var(cj) is obtained from the diagonal terms of [A]. 

In order to solve equation (6.44), the derivatives {R} and {F} must 

c i c i 

be found. Considering {F} , for example, the transformation is achieved by 

c i 

using the chain rule of differentiation, 


{ F} 


c i 


n _ 5b, 

E {Fk 7 ^ 
k=l b k *7 


(6.47) 


Since the random variable vector {b} is given by 

{b} = [*]{c} 

we have 





Hence, equation (6.47) becomes 
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(6.48) 



n 

E 

k=l 



ki 


Thus, prior to computing {F} from equation (6.48), all n vectors {F}, must 

c i D k 

be assembled and stored. 

It has been found that this diagonalization technique is highly beneficial 
when nonlinear computations are performed, since the eigenproblem in (6.40) 
can be solved once at the beginning and the result can be used for each load 
step. Thus, the double summations are reduced to single summations for each 
load step, without having to resolve the eigenproblem each time. 


6.7.2 Adjoint Method 

A distinct advantage in this first-order probabilistic finite element 
method is its ability to preselect certain displacement and stress responses of 
interest and thereby save computational effort by not calculating covariance 
response of others. Since the assumption has been made here that the random 
variables are independent from layer to layer, then the perturbation equations 
(6.43)— (6.46) must be solved for each layer as well. This could become 
expensive as the number of layers become large, therefore it would be nice to 
separate the computations from a dependency on the number of layers as much 
as possible. 

Liu et al. [34,40] introduced an adjoint method which, when used 
properly, partially achieves this goal. The method is used to calculate the 
displacement derivatives of the kth component of the generalized displacements 
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{A}. Using the chain rule of differentiation on {A}, we obtain 


where 




{D k } = <9A k /<9{A} 


(6.49) 


We can calculate {A} from Eq. (6.44) 

c i 


{&} = IkY'WL - <F} C ) = [kVw (6.50) 

c i c i c i c i 


Substituting (6.50) into (6.49) (note that the explicit derivative term ^ — is 


zero), we obtain 


aJ; = {dYikYYOc 

c i L i 


(6.51) 


The adjoint problem is defined by 


(K T ]{A k } = {D k } 


(6.52) 


Substituting (6.52) into (6.51), we arrive at the result 


A k = [K T HA k }[K T r 1 {f} c 

c i i 


or, finally, we obtain 
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(6.53) 


= (A k } T {f} c 

c i c i 

The result in equation (6.53) expresses the desired displacement derivative 
component in terms of the known force derivatives {f) c and the solution to 

k T 

the adjoint problem {A } . In solving the adjoint problem in equation (6.52), 
we note that the right hand side {D } simply becomes a Boolean vector, with 
unit value at the kth component. With this method, we simply select m 
displacement degrees of freedom (k = l,...m) of interest (usually only those 
necessary to calculate the stresses of interest), and solve equation (6.52) m 
times. Once m solutions are available, then equation (6.53) can be used to 
determine all ply level derivatives of random variables. Since the number of 
random variables is a function of the number of layers, equation (6.53) gives the 
derivatives for each ply without having to solve a structural level problem as in 
(6.44) for each set of random variables in each ply. Thus, we have succeeded in 
separating the solution of equation (6.44) from dependency on the number of 
layers in the model. Moreover, if the need to solve for a preselected m degrees 
of freedom exists, and m is much smaller than the number of nodal random 
variables times the number of layers, then substantial savings can result. One 
can easily see that for a large number of layers in the composite the adjoint 
method is very beneficial. 
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6.7.3 Assembly Savings Technique 

After gaining experience with the probabilistic methods in the shell 

program, it was discovered that considerable computational time is used in 

assembling the necessary residual force derivatives {F}, in equation (6.11). 

D i 

This expense is due to the fact that {F} b must be assembled independently for 

each "i" node in the model and for each of the P layers, or i*P assemblies for 
the solutions of equation (6.11). Recalling equation (5.41) for the residual force 
vector {F}, and substituting (5.42) we can express the integration of {F} b as 


{F} 



3 3 

E E 

r=l s=l 


[SDl] T {N} b + [SD2] T {M} b ° | J | w * w 
i i s r 


For the example of elastic materials, {N} b and {M} b become 


(6.54) 


< N >b. = [Al] b {SI} + |A2] b {S2} 


{M} b = [A2], {SI} + [A3]. {S2} 
i i i 


(6.55) 


where the general Green-Lagrange strains have been expressed as 

{<} = {SI} + ({S2} + C 2 {S3} (6.56) 

and terms with higher degrees than ( have been ignored. Equations (6.24) 
contain the derivatives [Al] b , [A2] b , [ A3] b for the case in which the ply 

thickness is a random variable. Each of these equations are in the form 
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(6,57) 


[A) b = PHOty 


where [D(£)] is some general function of ( and is the interpolation function 


n 


used in the expansion b(x) = E ip-(x)b-. In symbolic form the integration of 

i=l 1 1 

{F}^ involves terms of the form 


3 3 

{F} b = £ E 
D i r=l s=l 


[sdij a P(OK^){si} 


(6.58) 

s r 's 


for each random variable bj at the ith node. Applying the chain rule of 
differentiation to {F}^ , we obtain 


< F >b, = < F > b Hr = <fVi < 6 - 59 ) 


The result in (6.59) allows us to express (6.58) as 


= 


3 3 

E E 
r=l s=l 


[SDl] i [D(C)]{Sl} 


| J| w = {F} b ^j (6.60) 
s r 's 


In summary, the new approach is to assemble {F) b once for each layer, 
independent of the nodal random variable quantities b^, and then determine 
{F}^ after assembly for each node i by the use of equation (6.60). Even though 

{F}^ must still be assembled for each layer, the separate assemblies for each 
node are now successfully eliminated with considerable savings. 
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7. APPLICATIONS 


7.1 Introduction 

A number of problems have been solved to demonstrate various aspects 
and capabilities of the formulation developed in the present study. Verification 
of the accuracy of the first-order second— moment probabilistic finite element 
method (FOSM— PFEM) for all linear and nonlinear problems is carried out by 
comparison with Monte Carlo simulation results (for a description of the Monte 
Carlo method see Appendix B). Examples of linear problems are presented first 
to demonstrate modeling refinement requirements, sensitivities of the solution 
to various random variables, and effects of different levels of spatial correlation 
and input variance. Graphite— epoxy and metal matrix composite properties are 
compared as well. Geometric nonlinear examples involving post buckling of 
plates and shells are presented next, with sensitivity to individual random 
variables demonstrated throughout the load range. Finally, combined 
elastoplastic and geometric nonlinear examples are given for a tension specimen 
with a hole for both ARALL and Boron/ Aluminum composites. Experimental 
comparisons are given whenever possible. 

7.2 Linear Analyses 

The problem chosen as a comparison base for all examples in this section 
is a shallow spherical shell with simply supported boundary conditions and 
uniform external pressure. The problem description is given in Figure 7.1, and 
all input material properties, variances and levels of spatial correlation are 
given in Table 7.1. The coefficient of variation (COV) is defined as the ratio of 
the standard deviation to the mean for a given random variable. All examples 


99 



Two layer (0/90) 

R = 1000 in a = 50 in h = 1 in 


BC’s (Simply Supported) 

u = w =#2 = 0atx = a/2 

v = w=# 1 =0aty = a/2 

v=^ = 0aty = 0 u=# 1 = 0atx = 0 


Figure 7.1 Spherical shell under external pressure; out— of-plane 
displacement w measured at the center; stress measured at 

gauss point nearest center in 0 degree ply. 
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Table 7.1 


Material Properties and Statistics for Graphite— Epoxy 
Spherical Shell Problem 


Random 

Variable 

Mean 

Standard 

Deviation 

Coefficient of 
Variation 

V 

V 

E n 

15.75xl0 6 

7.8750xl0 5 

0.05 

25 

15 

E 22 

0.9091xl0 6 

4.5454xl0 4 

0.05 

15 

25 

°12 

0.4475xl0 6 

2.2376xl0 4 

0.05 

20 

20 

"12 

0.2223 

1.1116xl0 — 2 

0.05 

20 

20 

G 13 

0.4475xl0 6 

2.2376xl0 4 

0.05 

20 

20 

G 23 

0.3497xl0 6 

1.7487xl0 4 

0.05 

20 

20 

E fll 

31.0xl0 6 

1.55xl0 6 

0.05 

25 

15 

E f22 

2.0xl0 6 

l.OxlO 5 

0.05 

15 

25 

G fl2 

2.0xl0 6 

l.OxlO 5 

0.05 

20 

20 

"fl2 

0.2 

l.OxlO -2 

0.05 

20 

20 

"f23 

0.25 

1.25xl0 — 2 

0.05 

20 

20 

E 

m 

0.5xl0 6 

2.5xl0 4 

0.05 

25 

25 

"m 

0.25 

1.25xl0 — 2 

0.05 

25 

25 

FVR 

0.5 

2.5xl0 — 2 

0.05 

20 

20 

e 

0° ,90° 

2° 

— 

25 

25 

* 

6 

0.5 

.025 

0.05 

25 

25 


* Vindicates fiber orientation angle and 6 indicates ply thickness. 

The subscripts f and m stand for fiber and matrix; absence of a letter 
subscript indicates a ply— level property. The units are psi and inches 
where appropriate. 
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in this study used a COV of .05 except for the case of ply orientation angle, 
where the standard deviation was chosen to be 2 degrees. When spatial 
correlation was used, typical correlation lengths were around one half of the 
domain of the problem for each principal material direction. Only one quadrant 
of the shell was modeled in the interest of computational savings, even though 
the spatial correlation assumptions may not be symmetric across the shell. Of 
course, this is not a limitation of the procedures or computational model 
developed herein. 

7.2.1 Graphite— Epoxy Composite Shell 

The first example in this section involved the base line shell problem 
modeled with graphite— epoxy ply— level properties. The intent here is to 
illustrate both the agreement of the first— order second— moment probabilistic 
finite element method results with the Monte Carlo results, and to show the 
mesh refinement requirements in relation to the random field. Figures 7.2 
through 7.5 contain plots of out-of-plane displacement and stress for both 
mean and variance along the x-axis. A 2x2 mesh of nine— node Lagrange 
elements was used here with fully reduced (2x2) integration. The Monte Carlo 
solution is one that fully converged at 1500 simulations. While the agreement 
for the mean plots is good, the variance results are less satisfactory. Refining 
the mesh to a 4x4 nine-node element mesh, quite good results for the variances 
as well as the mean values were obtained (see Figures 7.6 through 7.9). This 
can be attributed to the fact that even though the 2x2 mesh was a reasonable 
one to discretize the deterministic equilibrium field equations, it was not refined 
enough to discretize the chosen random field. Thus the analyst must be 
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MEAN OF If DISPLACEMENT 



Figure 7.2 Mean center diplacement w along x— axis of spherical shell 
using all ply— level random variables and a 2x2 nine— node 
element mesh. 
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VARIANCE OF If DISPLACEMENT (xl.E-4) 



Figure 7.3 Variance of center displacement w along x— axis of spherical 
shell using all ply— level random variables and a 2x2 nine— node 
element mesh. 
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Figure 7.4 Mean o ^ stress along x-axis of spherical shell using all 
ply-level random variables and a 2x2 nine-node element mesh. 
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Figure 7.5 Variance of a stress along x— axis of spherical shell using all 
ply-level random variables and a 2x2 nine— node element mesh. 
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MEAN OF W DISPLACEMENT 



DISTANCE ALONG X-AXIS 


Figure 7.6 Mean center displacement w along x— axis of spherical shell 
using all ply— level random variables and a 4x4 nine— node 
element mesh. 
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VARIANCE OF If DISPLACEMENT (xl.E-4) 
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Figure 7.7 Variance of center displacement w along x— axis of spherical 
shell using all ply— level random variables and a 4x4 nine— node 
element mesh. 
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DISTANCE ALONG X-AXIS 


Figure 7.8 Mean stress along x— axis of spherical shell using all 
ply-level random variables and a 4x4 nine-node element mesh. 
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Figure 7.9 Variance of a stress along x— axis of spherical shell using all 
ply— level random variables and a 4x4 nine— node element mesh. 
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sensitive to both of these requirements when deciding whether a mesh is 
sufficiently refined. Generally, as the correlation lengths are increased, the 
effect on mesh requirements is that a coarser mesh can be used. For the rest of 
the spherical shell examples, this 4x4 nine-node mesh was selected. 

In Chapter 6 it was indicated that as the correlation lengths become 
larger, eventually the random field becomes fully correlated and equivalent to 
assuming a single random variable for that layer. This was called the "uniform 
variance" assumption, and the solution technique that results from this 
assumption involves considerably less computational expense. It is desirable to 
know when this method can be used versus the random field approach. Figures 
7.10 and 7.11 contain plots of displacement and stress variance versus the A^, 
correlation length normalized by the x— direction dimension of the shell. The 
X correlation length was selected to be five times the y— dimension of the shell 
and was held constant. The horizontal line in each figure is the uniform 
variance solution. From these results it is apparent that once both correlation 
lengths A^, and A^, are greater than four times their respective domain 
dimensions, the uniform variance solution can be used with equivalent results 
and considerably less expense. 

In all examples presented so far, a COV of .05 was used. A legitimate 
question remains as to how large the input variance (COV) can become before 
this first-order probabilistic method becomes inaccurate. Figures 7.12 and 7.13 
contain plots of the percentage difference of the first— order second— moment and 
Monte Carlo solutions versus the standard deviation for the ply angle and ply 
modulus random variables, respectively. From these results it can be 
concluded that for this problem if the ply angle standard deviation is less than 5 
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VARIANCE OF W DISPLACEMENT (xl.E-4) 
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Variance of center displacement w versus correlation length for 
spherical shell problem; comparison of uniform variance and 
random field results. 
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Figure 7.11 Variance of a stress versus correlation length for spherical 

shell problem; comparison of uniform variance and random 
field results. 
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% DIFFERENCE 



PLY ANGLE STANDARD DEVIATION (DEGREES) 


Figure 7.12 Percent difference between first-order second— moment method 
and Monte Carlo method for various standard deviations of the 
fiber orientation angle random variable. 
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COEFFICIENT OF VARIATION OF Ell 


Figure 7.13 Percent difference between first-order second— moment method 
and Monte Carlo method for various COV levels of 
random variable. 
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degrees, and the COV is less than .15, the percentage difference in the 
first— order second— moment and Monte Carlo solutions is less than 5%. This is 
reasonable, considering the fact that a 5 degree standard deviation for ply angle 
or a COV of .15 for is quite large for most composite (or isotropic for that 
matter) material variations. This is in agreement with Ang [23] and Liu et al. 
[33 — 40] , who stated that accuracy is maintained for a COV of .10 or less. Note 
that the first-order mean values here are nothing but the deterministic values, 
and that the agreement here is good without including the second— order 
perturbation effect. 

One very important benefit of the probabilistic finite element method is 
the ability to quantify the variations in the structural response caused by 
individual random variables. In the present study these random variables can 
include ply— level material stiffnesses or micro— level material stiffnesses, the 
latter evaluated with the aid of the Aboudi micromechanics model [64]. Figures 
7.14 and 7.15 illustrate both the combined and individual variances for 
displacement and stress response. It is interesting to note that for this 
particular shell problem, the w— displacement response is most affected by the 
ply thickness, ply angle, and E^ variables, while the stress in the 0 degree 
layer) is primarily influenced by the ply angle, with E^ and ply thickness 
variables much less significant. In Figures 7.16 and 7.17, micromechanics— level 
random variables were chosen. For the w— displacement variance, once again 
ply thickness and ply angle were important, along with FVR and E^ 
micro-variables. As for the stress, ply angle is still very dominant with E^,, 
FVR, and ply thickness secondary. It has been found that for these shell 
problems with graphite— epoxy materials, generally the dominant random 
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VARIANCE OF W DISPLACEMENT (xl.E-4) 


3.0 



Figure 7.14 Variance of center displacement w along x— axis of spherical 
shell showing the combined and individual effects of the 
ply— level graphite— epoxy random variables. 
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DISTANCE ALONG X-AXIS 


Figure 7.15 Variance of a stress along x— axis of spherical shell showing 

the combined and individual effects of the ply— level 
graphite— epoxy random variables. 
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VARIANCE OF W DISPLACEMENT (xl.E-4) 


3.5 



Figure 7.16 Variance of center displacement w along x— axis of spherical 
shell showing the combined and individual effects of the 
micro— level graphite-epoxy random variables. 
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DISTANCE ALONG X-AXIS 


Figure 7.17 Variance of a ^ stress along x-axis of spherical shell showing 

the combined and individual effects of the micro-level 
graphite— epoxy random variables. 
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variable is the ply angle, as would be expected due to the low stiffness of the 
matrix in comparison with the fibers. 

7.2.2 Metal Matrix Composite Shell 

To illustrate the differences in response when using metal matrix 
composite properties, the same shell was modeled using Silicon Carbide fibers in 
a Titanium Aluminide matrix. The input properties and statistical parameters 
used are listed in Table 7.2. Figures 7.18 and 7.19 illustrate the combined and 
individual variances for displacement and stress response with micro— level 
random variables. Now with the metal matrix properties, the w— displacement 
response is most affected by the matrix modulus E m and fiber volume ratio 
(FVR). This may be due to the fact that the mean FVR was selected to be .35, 
so if a higher mean FVR is input, fiber properties may dominate. As for the 
a stress, the longitudinal fiber modulus and FVR were the most 

important. Figures 7.20 and 7.21 contain the corresponding results but using 
ply— level random variables. Similar to the previous results, the E 22 and E^ 
moduli dominate the variance for the w— displacement response, and for the cr^ 
stress E^i is by far the most influential. 

7.3 Geometric Nonlinear Analyses 

Two problem types are discussed in this section. First, the 
graphite— epoxy spherical shell used in the previous section is analyzed deep into 
the postbuckling range, and second, the postbuckling of a flat non— stiffened 
graphite-epoxy panel with uniaxial compression is investigated. Experimental 
comparisons are made for the second example. 
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Table 7.2 


Material Properties and Statistics for Silicon Carbide 
Titanium Aluminide Spherical Shell Problem 


Random 

Variable 

Mean 

Standard 

Deviation 

Coefficient of 
Variation 

V 

V 

E n 

25.783xl0 6 

1.2892xl0 6 

0.05 

25 

15 

E 22 

18.683xl0 6 

9.3415xl0 5 

0.05 

15 

25 

G 12 

7.028xl0 6 

3.5140xl0 5 

0.05 

25 

25 

^12 

0.2706 

1.3528xl0 — 2 

0.05 

25 

25 

G 13 

7.028xl0 6 

3.5140xl0 5 

0.05 

25 

25 

G 23 

6.955xl0 6 

3.4774xl0 5 

0.05 

25 

25 

E fll 

50.7xl0 6 

2.535xl0 6 

0.05 

25 

15 

E f22 

50.7xl0 6 

2.535xl0 6 

0.05 

15 

25 

G fl 2 

21.3xl0 6 

1.065xl0 6 

0.05 

25 

25 

. 

S3 1 

0.19 

9.5xl0 — 3 

0.05 

25 

25 

^f23 

0.19 

9.5x10 3 

0.05 

25 

25 

E m 

12.3xl0 6 

6.15xl0 5 

0.05 

25 

25 

V 

m 

0.32 

1.6xl0 — 2 

0.05 

25 

25 

FVR 

0.35 

1.75xl0 — 2 

0.05 

25 

25 

* 

9 

0° ,90' 

2' 

— 

25 

25 

* 

6 

0.5 

2.5xl0 — 2 

0.05 

25 

25 


* 6 indicates fiber orientation angle and 6 indicates ply thickness. 

The subscripts f and m stand for fiber and matrix; absence of a letter 
subscript indicates a ply— level property. The units are psi and inches 
where appropriate. 
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Figure 7.18 Variance of center displacement w along x— axis of spherical 
shell showing the combined and individual effects of the 
micro— level Silicon Carbide Titanium Aluminide random 
variables. 
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DISTANCE ALONG X-AXIS 


Figure 7.19 Variance of a ^ stress along x— axis of spherical shell showing 

the combined and individual effects of the micro— level Silicon 
Carbide Titanium Aluminide random variables. 
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VARIANCE OF If DISPLACEMENT (xl.E-6) 



DISTANCE ALONG X-AXIS 


Figure 7.20 Variance of center displacement w along x-axis of spherical 
shell showing the combined and individual effects of the 
ply— level Silicon Carbide Titanium Aluminide random 
variables. 
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VABIANCE OF SIGMAX STRESS (xl.E3) 



DISTANCE ALONG X-AXIS 


Figure 7.21 Variance of a stress along x— axis of spherical shell showing 

the combined and individual effects of the ply— level Silicon 
Carbide Titanium Aluminide random variables. 
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7.3.1 Postbuckling of Spherical Shell 

All input properties and statistical parameters are the same as in section 

7.2.1 for the graphite— epoxy spherical shell. The difference here is that a 
geometric nonlinear solution is included, one that follows the nonlinear path 
past the snap through point (or limit point) and into the postbuckling range. 
The modified— Riks method [86,91] is used, as the regular Newton— Raphson 
method with load step control cannot trace this type of postbuckling curve. 

Once again, the Monte Carlo method is utilized to verify the first-order 
second— moment probabilistic results. Since the Monte Carlo method solves the 
nonlinear problem completely for each sample (or simulation) of the random 
variables, if the modified— Riks method is used in conjunction with the Monte 
Carlo method, each sample would result in a different set of load step sizes due 
to the self adjusting mechanism of the Riks method. The statistics used to 
estimate the mean and variance from the simulation results at each load step 
rely on all the variable responses residing at the same load value. For this 
reason the modified— Riks method was not used with Monte Carlo to verify the 
geometric nonlinear results. Instead, the Newton— Raphson method with 
constant load step size was used during this verification stage. Since the 
Newton— Raphson method could not pass the first limit point (zero— slope), then 
the solution was stopped there. Figures 7.22 through 7.25 contain the 
comparisons of the Monte Carlo and first— order second— moment solutions for 
mean and variance displacement and stress responses, and the agreement is very 
good. 
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UNIFORM PRESSURE 



MEAN OF CENTER DEFLECTION 


Figure 7.22 Mean center displacement w versus load for spherical shell 
using all ply— level graphite— epoxy random variables. 
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UNIFORM PRESSURE 



VARIANCE OF CENTER DEFLECTION (xl.E-4) 


Figure 7.23 Variance of center displacement w versus load for spherical 
shell using all ply-level graphite— epoxy random variables. 
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UNIFORM PRESSURE 



Figure 7.24 Mean of a stress versus load for spherical shell using all 
ply— level graphite— epoxy random variables. 
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UNIFORM PRESSURE 



VARIANCE OF SIGMAX STRESS (xl.E4) 


Figure 7.25 Variance of a stress versus load for spherical shell using all 
ply— level graphite— epoxy random variables. 
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Next, using the first-order second— moment method combined with the 
modified— Riks technique, the responses throughout the postbuckling range, 
including the limit points, are calculated and presented in Figures 7.26 and 
7.27. Figures 7.26a and b exhibit the w-displacement mean and COV. It is 
evident that at the limit points the COV is very large, almost 0.5. Figures 
7.27a and b contain the a ^ stress mean and COV. Again near a zero slope 
point on the stress curve the COV was quite large, almost 0.5. Figures 7.28 and 
7.29 illustrate the results in a different format. Here the squares indicate mean 
response while the stars are the mean plus or minus one standard deviation (one 
sigma) points. The influence of the limit points is more apparent in this 
format. 

The increase in variance at the limit points occurs since the buckling 
behavior of the structure is more sensitive to any changes in stiffness or load at 
these points. It should be noted that while the displacement COV begins at 
about .05 before the limit points, it becomes quite small after the limit points, 
settling to a value of about .01. As for the stress, the initial COV is .10, and 
after the limit points is still quite large, in the range of .07 to 1.6. 

7.3.2 Postbuckling of Flat Panel Under Axial Compression 

The problem under consideration here is a flat, rectangular 
graphite— epoxy panel loaded in axial compression. An experimental study was 
performed on a series of these panels by Starnes and Rouse [92]. Figure 7.30 
shows a typical panel with fixture and the resulting failure mode. The loaded 
ends of the panels were clamped by fixtures and the unloaded edges were simply 
supported by knife-edge restraints to prevent the panels from buckling as wide 
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Figure 7.26 Displacement postbuckling response of spherical shell using all 
ply— level graphite— epoxy random variables. 
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Figure 7.27 Stress postbuckling response of spherical shell using all 
ply— level graphite— epoxy random variables. 
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UNIFORM PRESSURE 


-4.5 



CENTER DEFLECTION 


Figure 7.28 Mean center displacement w and plus or minus one standard 
deviation points versus load throughout postbuckling range of 
spherical shell. 
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SIGMAX STRESS 


Figure 7.29 


Mean a stress and plus or minus one standard deviation 

points versus load throughout postbuckling range of spherical 
shell. 
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(a) Panel fixture 



(b) Panel failure mode 


Figure 7.30 Flat rectangular graphite— epoxy panel under axial compression 

[92]- 
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columns. The geometry, boundary conditions, and layup are given in Figure 
7.31, and the material properties and statistics are supplied in Table 7.3. 

A previous deterministic analytical study was performed on the panel 
(denoted as Panel C4 in [92]) by this author in reference [93]. Comparisons 
were made between analytical and experimental results. In general the 
comparisons were very good, even deep in the postbuckling range. In order to 
pass the critical buckling load, a geometric imperfection of a small percentage of 
the plate thickness (typically 1 to 5) times the normalized linear buckling mode 
was added to the original geometry of the panel. The purpose of the previous 
analysis was to study the effect of shear deformation on postbuckling response 
and failure prediction. The purpose of this analysis is to study the variability of 
the panel results. 

The probabilistic analysis of this panel assumed a fully correlated 
random field for each random function in each layer thus allowing the uniform 
variance solution to be used. An attempt was made to employ the random field 
techniques; however, the computational expense was too high due to eight 
random functions in 24 layers, 1625 degrees of freedom in the model, and a 
nonlinear analysis. The random field method took 8.3 hours per load step (on a 
Convex computer), whereas the uniform variance method only 5.8 minutes per 
load step. Since 13 load steps were required to reach the failure load and 
immense storage is required for the random field method with this many layers, 
it was concluded that for this size problem the uniform variance method was 
much more realistic. 

The end shortening postbuckling response is shown in Figure 7.32. The 
load P is normalized by the analytical buckling load P cr , and the end 
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Mesh Definition: 

72 nine-node Lagrange elements, uniform spacing with 12 
along x direction and 6 along y; 325 nodes and 1625 DOF. 


Figure 7.31 Geometry, boundary conditions, and layup for graphite-epoxy 
panel with axial compression. 
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Table 7.3 


Material Properties for Graphite— Epoxy Flat 
Panel Under Axial Compression Problem 


Random 

Variable 

Mean 

Standard 

Deviation 

Coefficient of 
Variation 

E n 

19.0xl0 6 

9.5xl0 5 

0.05 

E 22 

1.89xl0 6 

9.45xl0 4 

0.05 

G 12 

0.93xl0 6 

4.65xl0 4 

0.05 

^12 

0.38 

1.9xl0 — 2 

0.05 

G 13 

0.93xl0 6 

4.65xl0 4 

0.05 

o 

to 

CO 

0.25xl0 6 

1.25xl0 4 

0.05 

* 




9 

±45° ,0° ,90° 

2* 

— 

* 

6 

5.5125xl0 — 3 

4.167xl0 — 3 

0.05 


* 6 indicates fiber orientation angle and 6 indicates ply thickness. 
The units are in psi and inches where appropriate. 
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Figure 7.32 End— shortening postbuckling response of graphite— epoxy 

panel. 
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shortening deflection u by the analytical end shortening u cr at buckling. A 1% 
plate thickness geometric imperfection was used. The analytical results 
compare favorably with the experimental results. In addition, the plus or minus 
one standard deviation points indicate the variation in the data. It should be 
noted that in [93] , it was shown that one of the reasons for the good agreement 
here is the inclusion of the shear deformation in the element formulation. 
Figure 7.33 contains the out— of— plane deflection w near a point of maximum 
deflection normalized by the panel thickness t. Figure 7.34 shows the 
longitudinal surface strains e near a point of maximum out-of— plane deflection 
normalized by the analytical buckling strain e cr - In this comparison only the 
reduced integration Gauss point closest to the experimental strain gage was 
used to calculate the strains. Interpolation has been shown to improve 
agreement with the experimental results. For all three of these plots the 
COV is typically around 2%, except for the w displacement prior to and at 
buckling which was large. 

In order to understand the failure mode, it is necessary to see the 
nonlinear buckling mode shape. An analytical contour plot of the out-of-plane 
deflection at an applied load of 2.1 P cr is shown in Figure 7.35a. A moire fringe 
pattern photograph from reference [92] of the out— of— plane deflections at the 
same load is shown in Figure 7.35b. Both patterns indicate two longitudinal 
half-waves with a buckling— mode nodal line at panel midlength. 

In References 92 and 93, it was determined that the failure mode was 
primarily due to transverse shear stress (r xz in O' ply) along the midlength 
(nodal line) of the panel. In later work it has been found that although is 


142 


3 



0 ■) — i — i — i — i — i — i — i — i — i — | — i — i — i — i — i — i — i — i — i — | — i — i — i — i — i — i — i — i — i — 

0 12 3 


W/T 


Figure 7.33 Out-of-plane displacement 
graphite— epoxy panel. 
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Figure 7.34 Longitudinal surface strains on opposite outer surfaces of the 
graphite— epoxy panel near a point of maximum out-of— plane 
displacement. 
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(a) Contour plot of analytical 
results 


(b) Photograph of moire— fringe 
pattern 


Figure 7.35 Comparison of experimental and analytical out— of— plane 
displacement patterns at applied load of 2.1 P cr for the 
graphite— epoxy panel. 
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the primary failure mechanism, other stresses such as and also cause 
failure along the nodal line region in other layers. 

Proceeding with the analysis, Figure 7.36a shows a contour plot of the 
stress in the third layer of the laminate (a 0° ply) at an applied load of 2.1 
P cr - High compressive axial stresses occur along the longitudinal edges of the 
panel. The redistribution of the axial (a^) stresses for this 0° ply along the 
panel midlength is shown in Figure 7.36b for three different load levels. The y 
coordinate is measured from one side of the panel and normalized by the panel 
width b. The longitudinal membrane strain is redistributed to the edges of the 
panel after buckling. Typical COV values are .05 at P/P cr = 1.0, and range 
from .08 to .16 for the higher loads. The plus or minus one standard deviation 
points are shown in the figure. The material allowables for axial stress are 203 
ksi in tension, and 165 ksi in compression, so the axial stress at this location is 
well below this. 

The distribution of the transverse shear stress in the third layer of 
the laminate (a 0° ply) is shown in Figure 7.37a for a load of 2.1 P cr _ It is 
observed that high transverse shear stress develops along the nodal line of the 
panel. Figure 7.37b shows the redistribution of the stress for three different 
applied loads. At the experimental failure load of 2.1 P cr , the stress 
approaches the material allowable value of 9 ksi. The COV for the stress 
was typically .055 except at the critical buckling load when it reached a value of 
.12. 

It is of interest to determine the most significant random variables in the 
variance for the stress. Figure 7.38 contains the peak mean r ^ stress, 
along with the COV for the combined and individual random variables for 
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(a) Contour plot of axial 
stress distribution 


(b) Stress distributions across panel 
midlength 


Figure 7.36 Normal stress a distribution in the third-layer (0° ply) from 
the surface of the graphite— epoxy panel. 
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(a) Contour plot of transverse 
shear stress distribution 


(b) Stress distributions across panel 
midlength 


Figure 7.37 Transverse shear stress distributions in the third layer (0° 
ply) from the surface of the graphite-epoxy panel. 
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TAUXZ (KSI) 



Figure 7.38 Mean and coefficient of variation (COV) of r xz stress versus 

load throughout postbuckling range showing the combined and 
individual effects of the ply— level random variables for the 
graphite— epoxy panel. 
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increasing load values. It is observed that the stress variance is most 
influenced by G 13’ with the ply thickness effect increasing as the buckling load 
is passed and bending effects become more important. 

Even though for this graphite-epoxy composite panel a first ply failure 
does not represent overall panel failure, reliability analysis was performed for 
first ply failure. For simplicity the maximum stress failure criterion was 
selected, although other failure criteria may be more suitable such as maximum 
strain or a quadratic polynomial criterion, such as the Tsai— Wu criterion. In 
this analysis, since the material allowables’ accuracy was a question, a very 
simple example reliability calculation was determined suitable in which the 
failure criterion was based on exceeding the transverse shear strength of 9 ksi. 
The strength and stress were selected to be Normal distributed, with the stress 
mean and COV from the finite element results and the strength COV assumed 
to be 0.10. Appendix B describes the reliability theory used here [23], Figure 
7.39 is a plot of the probability of safety curve (probability of non— failure of 
first ply) versus load for the failure mode. Based on the first— order 
second-moment probabilistic finite element results and the assumed strength 
statistics, approximately 100% reliability against first ply failure exists at a 
load step of 1.65 P cr . It is believed that the reliability analysis for first ply 
failure used in conjunction with a progressive failure analysis, in which the 
damage due to local failures is progressively accumulated, would be a very good 
criterion for design of these panels to exclude any failures. Another option 
would be to determine the system reliability, in which individual ply 
reliabilities are combined into the overall laminate structure reliability. 
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Figure 7.39 Reliability against first ply failure due to r xz stress versus 

load throughout postbuckling range of the graphite-epoxy 
panel. 
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7.4 Material Nonlinear Analyses 

Material nonlinearity has been included in a fashion such that combined 
material and geometric nonlinearity can be studied. The classical example of an 
isotropic cylindrical shell roof under self weight is analyzed and compared with 
results in the literature for all three cases of geometric nonlinearity, material 
nonlinearity, and combined nonlinearities. The same example is then used to 
validate the mean and variance statistics of the first— order second— moment 
probabilistic method by comparison to Monte Carlo results for the combined 
nonlinearity case. 

Two material nonlinear composites were then selected for study. First, 
an ARALL laminate (patented by Alcoa), composed of aramid epoxy layers in 
between layers of isotropic aluminum, was analyzed. Second, a single layer 
Boron/Aluminum composite is modeled using the orthotropic plasticity 
formulation. Both materials were then employed in a model of a tension 
specimen with a hole. 

7.4.1 Cylindrical Shell Roof Under Self Weight 

The cylindrical shell roof under self weight problem described in Figure 
7.40 is a classical test example in the literature. The shell has free longitudinal 
edges and is supported at both ends on rigid diaphragms. Ideal plasticity is 
assumed for the material nonlinear behavior. The material properties and 
statistics are given in Table 7.4. 

In order to validate the plasticity model shell formulation, the 
deterministic results for this problem are compared to those obtained by Ramm 
and Sattele [94]. In their work a 3x3 mesh of bicubic 16— node degenerated shell 
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Figure 7.40 Cylindrical shell roof under self weight problem description. 
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Table 7.4 


Material Properties and Statistics for Cylindrical Shell Roof Under 

Self Weight Problem 


Random 

Variable 

Mean 

Standard 

Deviation 

Coefficient of 
Variation 

V 

V 

E 

21000 

1050 

0.05 

2440 

3800 

V 

0.0 

— 

— 

— 

— 

* 

a Y 

4.2 

0.21 

0.05 

2440 

3800 


o 

Material is isotropic; v is not considered random. Units are in N/mm and 
mm where appropriate. 

* Oy stands for yield stress. 
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elements was used with a total of 441 degrees of freedom. The cubic element 
uses 4x4 Gauss integration in the plane of the element and a 7 point Simpson’s 
rule integration in the thickness direction. Figure 7.41 contains the results 
obtained by Ramm and those for the present study for the center node 
deflection w^ at the free edge. Various mesh refinements of biquadratic 9-node 
elements were used here to study agreement with the Ramm solution. Full 
integration for the biquadratic element consists of 3x3 Gauss points. This 
element is known to have locking (over— stiffening) problems when used to 
model thin shells, and this is typically remedied by using either fully reduced 
integration on all terms (2x2) or selective reduced integration in which 3x3 
integration is used for the in— plane terms and 2x2 integration for the transverse 
shear terms. Studying the figure, the 4x4 element mesh with full 3x3 
integration obviously exhibits locking for the combined nonlinear results and 
thus is too stiff. Using fully reduced integration, the geometrical nonlinear 
elastic, and linear elastoplastic solutions agree well with the reference. However 
the combined nonlinear solution is too soft. By refining the mesh the locking 
effects are reduced such that for a 6x6 biquadratic mesh with 3x3 integration 
the results compare very favorably with those obtained by Ramm. The degrees 
of freedom in the three biquadratic meshes are as follows: 4x4 mesh — 405 

d.o.f., 5x5 mesh — 605 d.o.f., and 6x6 mesh — 845 d.o.f. Eight gauss integration 
points through the thickness were used. In addition, it is of interest to note the 
stiffening effect due to the geometrical nonlinearity and the softening effect due 
to the elastoplasticity. 

The cylindrical shell problem is used next to validate the first-order 
second— moment probabilistic finite element mean and variance response by 
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Figure 7.41 Center node displacement w at location A of cylindrical shell 
roof versus load for various types of nonlinearity. 
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comparison with Monte Carlo results. The input statistics are given in Table 
7.4, and to save computation time the 4x4 model with 2x2 integration is used. 
The comparisons are given in Figures 7.42 — 7.45, and are made for the 
combined nonlinear solution, with 1500 Monte Carlo simulations. The 
deterministic solution yields at a load of 1.4x10 , and five load steps past yield 
are shown. More load steps are not used as the expense of the Monte Carlo 
solution prohibits this. The Monte Carlo solution required 43.4 hours CPU time 
while the first-order second— moment method only used 5.3 minutes. Obviously 
the perturbation method is computationally advantageous. 

7.4.2 ARALL Laminate Tension Specimen with Hole 

ARALL laminates are high strength hybrid composites for aerospace 
applications developed by Alcoa. Figure 7.46 illustrates the concept of bonding 
thin sheets of high strength aluminum alloys using high strength aramid fibers 
in a special epoxy resin. The benefits include significant increases in fatigue and 
fatigue crack growth properties over monolithic aluminum, and the outer 
aluminum layers provide impact damage and moisture protection that would be 
a problem for typical fiber composite materials. In addition, increases in 
strength and lower densities are achieved as compared to monolithic aluminum 
[96]- 

It is desired to model the nonlinear structural behavior of this hybrid 
composite. Figure 7.46 also shows the description of the analytical model used 
to represent the stiffness of this material. The aramid epoxy layers are divided 
into fiber— rich and resin— rich layers with stiffness given for each sublayer. 
Table 7.5 gives the properties and statistics for the aluminum, aramid epoxy 
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VERTICAL LOAD P (xl.E3) 



Figure 7.42 Mean center node displacement w at location A versus load for 
the cylindrical shell roof under self weight. 
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VERTICAL LOAD P (xl.E3) 



VARIANCE OF W DEFLECTION AT A 


Figure 7.43 Variance of center node displacement w at location A versus 
load for the cylindrical shell roof under self weight. 
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VERTICAL LOAD P (xl.E3) 



Figure 7.44 Mean a stress at location B versus load for the cylindrical 
shell roof under self weight. 
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VERTICAL LOAD P (xl.E3) 



Figure 7.45 Variance of stress at location B versus load for the 
cylindrical shell under self weight. 
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Figure 7.46 ARALL laminate layup and geometry. 
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Table 7.5 

Material Properties and Statistics for ARALL— 1 Laminate 

Constituents 


Random 

Variable 

Mean 

Standard 

Deviation 

Coefficient of 
Variation 

Aluminum (7075— T6L) 



E 

10.4x10° 

5.2xl0 5 

0.05 

V 

0.3 

1.5xl0 — 2 

0.05 

* 

°Y 

7.8xl0 4 

3.9xl0 3 

0.05 

* 

5 

1.2xl0 — 2 

6-OxlO -4 

0.05 

Aramid Edoxv fiber— rich layers 



E n 

12.549x10° 

6.2745xl0 5 

0.05 

E 22 

0.76525xl0 6 

3.82625xl0 4 

0.05 

G 12 

0.28955xl0 6 

1.44775xl0 4 

0.05 

v 12 

0.3458 

1.729xl0 — 2 

0.05 

G 13 

0.28955xl0 6 

1.44775xl0 4 

0.05 

G 23 

0.26462xl0 6 

1.3231xl0 4 

0.05 

* 

e 

0° ,90° 

2° 

— 

* 

5 

5.6xl0 — 3 

2.8x10^ 

0.05 

Aramid Edoxv resin— rich layers 



E 11 

2.1972x10° 

1.0986xl0 5 

0.05 

E 22 

0.48219xl0 6 

2.41095xl0 4 

0.05 

G 12 

0.15717xl0 6 

7.8585xl0 3 

0.05 

"l2 

0.3749 

1. 8745x1 0 -2 

0.05 

G 13 

0.15717xl0 6 

7.8585xl0 3 

0.05 

G 23 

0.15576x10^ 

7.7880xl0 3 

0.05 

* 

9 

0° ,90* 

2* 



* 

6 

1.416xl0 — 3 

7.08xl0 — 5 

0.05 


* <7y indicates yield stress, 0 indicates fiber orientation angle, and d indicates 
ply thickness. 

Units are in psi and inches where appropriate. 
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fiber— rich, and aramid epoxy resin— rich layers. Experimental tension test 
results [96] are compared to the analytical results in Figure 7.47. From the 
figure it is observed that the aramid epoxy behavior is linear and the analytical 
linear comparison is very good. The 7075— T6(L) aluminum behavior is elastic 
perfectly-plastic and the analytical model with a yield stress of 78 ksi agrees 
very well except near the point of first yield. As for the ARALL-1 results (-1 
indicates 7075-T6(L) aluminum is used) the analytical model with ideal 
plasticity for the aluminum layers and linear elastic aramid epoxy layers 
generally exhibits the same behavior as the experimental results except the 0 
degree laminate analytical model underpredicts the stiffness after yield and the 
90 degree laminate model overpredicts the stiffness after yield. For the purpose 
of this example the analytical model is considered acceptable and will be used to 
study the mean and variance response of an ARALL tension specimen with a 
hole. 

Figure 7.48 shows the finite element model and dimensions of the tension 
specimen with a hole problem. The same material properties and material 
model from the previous discussion are used in this problem. The probabilistic 
analysis assumed a fully correlated random field for each random function in 
each layer, thus allowing the uniform variance technique to be used here. Again 
the computational expense of discretized random fields in each layer in a 
nonlinear analysis forced this assumption. Figure 7.49 contains the mean and 
standard deviation of the longitudinal e strain at the hole edge (point A in 

yy 

Figure 7.48) for the case where all aramid epoxy layers are aligned at 90 degrees 
to the load. The figure also shows the breakdown of standard deviations for all 
the significant random variables. Since the fibers are 90 degrees to the load and 
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TENSILE STRESS (KSI) 



Figure 7.47 ARALL tension test experimental and analytical comparisons. 
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BC’s: at x = 0 u = 0 

at y = 0 v = 0 

at y = L u = 0 


Figure 7.48 Finite element model and dimensions of ARALL tension 
specimen with hole. 
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Figure 7.49 ARALL tension specimen with hole; mean and standard 
deviation of e strain at the hole edge (point A) versus load 

for the case where the fibers are 90° to the load; included are 
the combined and individual effects of the ply-level random 
variables. 
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to the strain e , then the aluminum properties tend to dominate. It is 

yy 

interesting to note that even though no bending occurs in this problem, the ply 
thickness of the aluminum layers is dominant after yield. The aluminum yield 
stress and elastic modulus are also important. Figure 7.50 contains similar 
results except now the fibers are aligned with the loading direction. While the 
aluminum yield stress and ply thickness random variables are still significant, 
the aramid and ply thickness random variables are now equally important. 
These results illustrate the role the individual random variables play in the 
total variability of this type of ARALL structure. 

7.4.3 Boron/ Aluminum Tension Specimen with Hole 

A Boron/ Aluminum laminate was selected to illustrate the use of the 
macroscopic orthotropic plasticity formulation. The same problem dimensions 
(except for thickness) were used as in the last example, however, as shown in 
Figure 7.51, a different mesh was used that placed gauss points along the 
x— axis. Rizzi, Leewood, Doyle, and Sun [75] conducted an experimental and 
analytical study of this specimen, and provided experimental measurements for 
the orthotropic elastic constants as well as the a^ values in the yield criterion 
and the hardening parameters in the isotropic work hardening model. These 
values are all stated in Table 7.6 and are used in the present analytical model. 
It should be noted that the a^ values used in this study differ from those given 
in the reference by a factor of 2/3 due to a minor difference in the formulations. 
Figure 7.52 contains a comparison of the analytical results from the present 
study and experimental results from reference [75] for the longitudinal strain 
e along a radial line (x-axis) 90 degrees to the loading. The agreement is 
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Figure 7.50 ARALL tension specimen with hole; mean and standard 
deviation of e ^ strain at the hole edge (point A) versus load 

for the case where the fibers were alligned with the load; 
included are the combined and individual effects of the 
ply— level random variables. 
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Figure 7.51 


P 

♦ ♦ ♦ 



Thickness = 0.0795 in 

BC’s: at x = 0 u = 0 

at y = 0 v = 0 

at y = L u = 0 


Finite element model and dimensions of Boron/ Aluminum 
tension specimen with hole. 
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Table 7.6 


Material Properties and Statistics for Boron/ Aluminum Laminate 


Random 

Variable 

Mean 

Standard 

Deviation 

Coefficient of 
Variation 

E n 

29.4xl0 6 

1.47x10® 

0.05 

E 22 

19.1xl0 6 

9.55xl0 5 

0.05 

^12 

7.49xl0 6 

3.745xl0 5 

0.05 

U 12 

0.169 

8.45xl0 -3 

0.05 

G 13 

7.49xl0 6 

3.745xl0 5 

0.05 

G 23 

7.49xl0 6 

3.745xl0 5 

0.05 

* 

Oy 

13.5xl0 3 

6.75xl0 2 

0.05 

_* 

H 

60.0xl0 3 

3.0xl0 3 

0.05 

♦ 

0 

0° 

2.0° 

— 

* 

6 

7.95xl0 — 2 

— 

0.05 


* ay indicates yield stress, H indicates hardening modulus, 6 indicates ply 
orientation angle, 6 indicates ply thickness. 

The values of the a- j constants in the yield criterion are: 

§ a n = 0.001 j a 22 = 1.0 § a 12 = -0 01 

3 _ 3 _ 3 . q 

2 a 44 ” 2 a 55 ~ 2 a 66 “ y 


a y \ 1 / \ 

The hardening model used was Y(a) = H[a + [— — ] ] ' 

H 


A = 5.8 

Units are in psi and inches where appropriate. 
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LONGITUDINAL STRAIN (xl.E-6) 



Figure 7.52 Analytical and experimental comparison of the longitudinal 
strain along a line 90° to the loading for the 

Boron/ Aluminum laminate tension specimen with hole. 
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slightly worse than that obtained in [75], but is probably due to the difference 
in element formulations and the classical incremental plastic stress routine used 
versus the radial return algorithm used here. Yielding occurs after 1000 lbs, 
and the agreement worsens as the loading is increased. However, the results are 
still considered quite good. 

Using the random variable statistics stated in Table 7.6, the first-order 
second— moment probabilistic method was used to evaluate the mean and 
variance of the e strain response. Once again the probabilistic analysis 

yy 

assumed a fully correlated random field for each random function which allowed 
the uniform variance technique to be used here. Figure 7.53 shows the 
analytical mean e strain for the 2500 lb and 1000 lb load values with the plus 

yy 

or minus one standard deviation points included. It is obvious that the 
sensitivity of to the random variables increases both with the load and as 
the location moves closer to the hole. Figure 7.54 is a plot of both the mean 
and standard deviation of the e strain at the location A on the model versus 
load. The breakdown for each random variable is presented as well. Since only 
a single layer is used, then the ply thickness could not be considered a variable 
here. The most significant random variable is the plastic hardening modulus H, 
with E 22 and the yield stress important as well. Note that E 22 is significant 
since the fibers are 90 degrees to both the loading direction and to e yy . This 
example could obviously be extended to include the a- • plastic yield coefficients 

J 

and the hardening parameter A as random variables since they are also 
experimentally measured quantities with uncertainties. 
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Figure 7.53 Analytical results for the longitudinal strain e along a line 

90° to the loading for the Boron/ Aluminum laminate tension 
specimen with hole. 
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Figure 7.54 Boron/ Aluminum tension specimen with hole; mean and 

standard deviation of e ■ longitudinal strain at the hole edge 

(point A) versus load; included are the combined and 
individual effects of the ply— level random variables. 
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8. CONCLUSIONS AND RECOMMENDATIONS 
8.1 Summary and Conclusions 

A probabilistic finite element analysis procedure for laminated composite 
shells is developed. Full geometric nonlinearity for large deformation and 
rotation and rate— independent anisotropic plasticity are included. A 
degenerated 3— D laminated composite shell element with first-order shear 
deformable kinematics and a total Lagrangian finite element formulation is used 
in the deterministic analysis. The first-order second-moment technique for 
probabilistic finite element analysis of random fields is employed to determine 
mean and variance of displacement, strain, and stress fields. Random variables 
built into the model include ply stiffnesses, orientation angles, and ply 
thicknesses. Fiber and matrix stiffnesses and volume ratios can be selected as 
random variables with the use of the Aboudi micromechanics model. Monte 
Carlo simulation was used to verify selected results. 

Many problems were investigated either to verify the second-moment 
method’s accuracy or to investigate and quantify variability in certain 
structures. It was concluded very early that the second— order perturbation of 
the second— moment method required too much computational expense and 
storage and returned only a slight correction to the mean values. By comparing 
the results with the Monte Carlo method, it is concluded that the first— order 
second— moment method for estimating structural response mean and variance is 
quite accurate as long as the input coefficient of variations are less than 0.15. 
The method maintains this accuracy when a significant number of independent 
random variables are assumed, as is typically the case in layered composite 
problems. A few large degree of freedom and large number of layer nonlinear 
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problems were studied to test the probabilistic method’s ability to deal with 
practical (aerospace industry) size models. The random field techniques for the 
larger models become far too costly both from computational and storage 
viewpoints. However, the assumption of a fully correlated field for each random 
variable in each layer (uniform variance) led to more realistic computational 
expense. 

The inclusion of transverse shear deformation proved to be critical in 
modeling laminated composites, especially into the postbuckling range. It was 
demonstrated that the modified Riks arc length method works quite well with 
the second— moment probabilistic method and allowed mean and variance 
calculations to be made beyond zero-slope limit points, which often exist in 
shell structures. 

As for material nonlinear problems, the radial return algorithm was 
installed in a manner such that combined geometric and material nonlinear 
problems can be solved quite efficiently. The plasticity analysis is performed 
here in combination with the geometric nonlinearity and resulted in very little 
increase in iterations per load step. ARALL and Boron/Aluminum plasticity 
problems were investigated and the variability of these composites were 
quantified for a tension specimen with a hole. 

An approximate reliability calculation against first ply failure was made 
for a composite panel loaded in axial compression into its postbuckled state. By 
assuming the stress and strength to have Normal distributions, the mean and 
variance of the stress could be used directly in a linear (maximum stress) 
performance function to estimate probability of safety. 
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8.2 Recommendations 

The natural and most important step to be added to the computational 
procedure developed herein is to efficiently and optimally integrate it with the 
first and second— order reliability estimation methods. Thus calculation of the 
mean and variance response would become optional if reliability computations 
are desired. The sensitivity derivatives already computed in the program would 
simply be used directly to compute the safety indices, and in this way the 
original random variables can be assumed to be any known distribution 
functions with specified mean, variance, and correlation. The loss in accuracy 
in calculating small probabilities of failure by using the mean and variance 
response directly would be avoided in this manner. 

With the improved reliability algorithm installed, more detailed and 
accurate laminated composite reliability calculations could be made. More 
generally accepted failure criteria, such as the Tsai— Wu, Tsai— Hill, and 
maximum strain criteria could be used. Laminate failure could be studied by 
combining the individual first ply failure probabilities into a system reliability 
problem. More detailed stress analysis may be required, leading to a 
global— local approach combining other theories such as the layer-wise theories 
of Reddy [97]. Stiffener elements could be added to model the effects of 
stiffened composite plates and shells. Also, design optimization could be 
performed by optimizing the reliability index, with selected random variables 
used as design variables. 

The generality of the formulation could be improved by broadening the 
range of loading and problem types. The addition of thermal loads, random 
loads and geometry, and transient analysis are obvious extensions. The 
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inclusion of plasticity and other nonlinear behavior at the constituent level of a 
micro mechanics theory is also a possibility. Due to the realization that the 
mechanics of the interphase region in a metal matrix composite is very critical, 
a micromechanics theory which includes this region is essential. 
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APPENDIX A 


In this appendix, a few relevant terms are defined and the basic method 
for determining reliability is described. This brief review is only meant as an 
aid to the reader. Most of the definitions involved in the initial part of this 
section are based on material from reference [56]. 

In general terms, a random variable is one whose value is uncertain or 
undetermined. For this reason problems without random variables are often 
considered to be deterministic. The distribution function controls the 
probability of a random variable having certain values. Two distribution 
functions are of major importance: the cumulative distribution function (CDF), 
and the probability density function (PDF). The CDF of a random variable U 
is given by 

Fjjfu) i P[U < u] (A.l) 

which means the probability that the random variable U is less than or equal to 
some deterministic value u. The PDF is defined as 


f u(“) 1 


du 


so that 


(A.2) 


P[U < u] = F^u) = /“ fyMdz 

— QD 


The weighted averages or moments of random variables can also be used 
to describe their distribution. The expected value of a function r(U) is defined 
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as 


E[r(U)] = j" r(u)f u (u)du 

— QD 


(A.3) 


Using this definition, important moments can be defined such as mean and 
variance. The first moment, called the mean, is the central measure of the PDF 
and is given by 


E[U] = Mu = | u %( u )du 


(A.4) 


The second moment (central), called the variance, is a measure of the dispersion 
or spread of the distribution from its mean and is defined as 

Var(U) i E[U - ^] 2 = |" (u - >* u ) 2 f u (u)du (A 5) 

— OD 


The standard deviation and coefficient of variation are parameters which are 
often utilized. The standard deviation is given by 

Gjj = N /Var(U) (A. 6) 

and the coefficient of variation (COV) by 
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(A.7) 


°\J 

COV = — 

Mu 


When two random variables have a probabilistic relationship, this 
becomes a family and the relationship is defined by their joint distribution 
functions. For example, the joint CDF of two random variables U and V is 
defined as 

F uv (u,v) = P[U < u,V < v] (A. 8) 


where the comma indicates the intersecting areas of the individual CDFs. 
Moments can also be defined for two random variables. The expectation of a 
function of U and V (first moment) is given by 


E[r(U,V)] = J J r(u,v)f uv (u,v)dv 

OD — QD 


(A.9) 


The second central moment, or covariance, 


Cov(U,V) E E[(U - % )(V - ^>1 = j" j" (u - /^)(» - /Xyjfuyfu.vjdv 

—a) “OD 

(A. 10) 


and corresponding correlation coefficient, 


191 


P UV E 


Co v ( U,V) 
a U a V 


are dimensional and dimensionless measures of linear dependence between the 
two random variables. 

A random function is an extension of these ideas in which this function 
varies with one or more variables but for specific values of the variables its 
value is uncertain. If the variable is restricted to time only, then the random 
function is called a random process. If the variables are only spatial 
coordinates, then the function is called a random field. In the present study 
random fields are the focus, however it is also shown that if a random field 
becomes very highly correlated, this reduces to or becomes equivalent to a 
single random variable. A random function has an infinite number of 
distribution functions, often referred to by "orders", which imply moments. For 
example, the second— order CDF of a random field U(x), where x refers to 
spatial coordinates, evaluated at x 1 and X2 is given by 


F U u ( u i> u 2 ,x l ,x 2^ = P [ U 1 = U ( x l) - u i ,U 2 = U ( x 2^ - u 2^ 

1 2 

(All) 

A homogeneous (or stationary if x refers to time) random function U(x) of order 
two is one whose second— order CDF and PDF are dependent only on the 
difference (x^— X2) and not on the actual values of x^ and X2. When the 
function has a constant mean and an autocovariance dependent only on this 
difference (x^— X2), then it is termed homogeneous in the wide sense. 
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Autocovariance refers to the covariance of a single random function evaluated 
at two points, given mathematically by 

rOD rCD 

Cov(U(x 1 ),U(x 2 )) = % u (u 1 ,u 2 ;x 1 ,x 2 )du 1 du 2 

J J 12 

— OD — OD 

(A.12) 

In the present study only random fields which are homogeneous in the wide 
sense have been considered. However, the methods can easily be applied to 
inhomogeneous fields in which the mean is a function of the actual spatial 
coordinates, for example. 

Having defined some basic probabilistic principles and terms, the next 
step is to describe the framework typically used to estimate reliability. The 
description to follow is modeled after the original ideas developed by Hasofer 
and Lind [78] and discussed by Wirsching and Wu [83]. Let g(lL) = 0 be the 
limit state function in which IL are the random variables. Each IT is 
transformed to a reduced coordinate u^ according to 

where (/z.,cr.) are the mean and standard deviation respectively of IL. Equation 
(A. 13) is then substituted into g(lL) so that the limit state function is now 
expressed in terms of the reduced coordinates, g^(u-). The generalized safety 
index /? is defined as the minimum distance from the origin of the reduced 
coordinates to the limit state surface. Thus in mathematical terms the problem 
becomes the following constrained minimization problem 
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P = min 



(A.14) 


subject to the constraint g^(uj) = 0. The design point is also defined as the 
point on g^( u j) = 0 closest to the origin. Typically, various optimization 
schemes are employed to solve the problem stated above. 

The probability of failure can be easily determined using the safety 
index (3 and the standardized Normal cumulative distribution tables (<p) as 


If g(U-) is linear in the U-, and all Uj are Normal, then Pj is exact. Otherwise 
Pf is only an approximation. 


and Fiessler [79], Chen and Lind [80], and Wu [81-85] which provide 
improvements to the estimate of P^. These improvements allow for other 
distributions than the Normal, and also nonlinear limit state functions. The 
recent Wu method has proven to be the most accurate and also the most 
complex. 

In the process of solving the constrained minimization problem, 
derivatives of the limit state function with respect to random variables are 
required. Since the limit state function is typically a function of the structural 
response (e.g. displacement, strain, stress, etc.), then by the chain rule of 
differentiation, derivatives of the response with respect to the random variables 
are needed. Since the latter derivatives are already determined in the process of 
calculating the variance using the second— moment probabilistic finite element 


p f = *(-&) 


(A.15) 


Extensions to the above basic method have been developed by Rackwitz 
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method, then they can now be used directly to estimate reliability. In this way 
an efficient procedure for reliability estimation incorporating the probabilistic 
finite element method is achieved. 
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APPENDIX B 


The Monte Carlo simulation method is considered to be a mature 
subject, with widespread use since the onset of rapid computers. The term 
"simulation", as used in the present context, is the technique of using a 
computer to evaluate a given deterministic model numerically. Monte Carlo 
simulation, in a general sense, can be defined as any simulation involving the 
use of random numbers for solving certain stochastic problems. In the following 
paragraphs, a brief description of the Monte Carlo simulation method is 
presented. Most of the content of this review is based on the material in 
reference [1]. 

Using the Monte Carlo simulation procedure, the computer is used to 
generate n independent statistical samples for each random variable, which are 
then fed into the model. Each sample can be thought of as an independent 
deterministic experiment, which is processed by the model to yield the results of 
the experiment. Each "sample" is drawn from a pre-selected probability 
distribution, so that the sample distribution is appropriate. After the 
simulation process, the output data is statistically analyzed to estimate the true 
characteristics of the model. 

Many probability distribution functions exist in the statistical literature. 
In this study the random variables were all assumed to follow the Normal 
distribution, whose probability density function (PDF) is defined as 



(B-l) 
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where and <Jjj are the mean and standard deviation of the random variable 
U. In the process of generating random numbers, the Uniform distribution 
becomes important. This distribution gives an equal probability of any number 
within the prescribed interval, and its PDF is defined as 

%(u) = 5-=“a if a < u < b (B.2) 

0 elsewhere 

Therefore once independent random numbers have been generated from the 
Uniform distribution, they can be transformed into random variables from any 
other distribution. Many random number generators and distribution sampling 
schemes exist in the literature, and the reader is therefore referred to [1]. 

When correlation between random variables exist, such as in the case of 
a random field, then the independent random variable samples produced from 
the random number generator/ sampling scheme must be transformed into 
correlated random variable samples. In this study a method proposed by 
Shinozuka [6] is used, in which the Choleski decomposition of the covariance 
matrix is utilized to perform this transformation. 

Since in the present work the mean and variance are the distribution 
parameters of the structural response under study, then it is necessary to 
estimate these moments using Monte Carlo techniques. With this in mind, the 
sample mean is defined as 
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(B.3) 


U(n) 


n 


E 

i=l 

n 


u. 

1 


such that U(n) is an unbiased estimator of the actual mean ^ . The sample 
variance is defined as 


E [ Uj - U(n)] 2 

s 2 (») * 1-1 n _ ! (B.4) 

2 2 
such that s (n) is an unbiased estimator of the actual variance cr^. The strong 

law of large numbers guarantees if a sufficiently large sample size n is taken, 

that U(n) 2 /Ay will be true. Thus it is obvious that the Monte Carlo 

simulation method requires a large sample size to be accurate, so that when this 

method is combined with the finite element method, considerable computational 

expense is the result. For this reason the Monte Carlo method is not attractive 

in performing probabilistic finite element analysis. However, it is considered 

quite accurate when large enough samples are taken, and for this reason it is 

used to check selected results of the second— moment probabilistic finite element 

method. 
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